///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/all_pair_jast.cpp
///
/// \brief   implementation file for the all pair jastrow coefficient object
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<vector>
#include<complex>
#include<string>
#include<algorithm>
#include<cmath>
#include<cassert>

#include<formic/fqmc/all_pair_jast.h>
#include<formic/wfn/configuration.h>
#include<formic/exception.h>
#include<formic/all_data.h>
#include<formic/command/command.h>
#include<formic/mpi/interface.h>
#include<formic/lapack/interface.h>
#include<formic/numeric/numeric.h>
#include<formic/output/output.h>
#include<formic/timing/timing.h>
#include<formic/fqmc/apj_fortran.h>
#include<formic/fqmc/jastrow_utility.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief builds the dense alpha-alpha or beta-beta jastrow factor array from the compact array
///
/// \param[in]      jaa      size _naa_var.  The compact array.
/// \param[out]     daa      size 4*_norb*_norb (simple = _norb*_norb).  On exit, the dense array.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::fill_dense_aa(const S * const jaa,
                                                       S * const daa) {

  if (_simple) {

    for (int i = 0; i < _norb; i++)
    for (int j = i; j < _norb; j++) {
      const S val = jaa[ formic::cmpd_pair_index_incl(i, j, _norb) ];
      daa[i+_norb*j] = val;
      daa[j+_norb*i] = val;
    }

  } else {

    std::fill(daa, daa + 4*_norb*_norb, formic::zero(S()));
    const S * ptr = jaa;
    for (int j = 0; j < _norb; j++)
    for (int i = 0; i < j; i++) {

      S * const pij = daa + 4*(i+_norb*j);
      pij[0] = ptr[0];
      pij[1] = ptr[1];
      pij[2] = ptr[2];
      pij[3] = ptr[3];

      S * const pji = daa + 4*(j+_norb*i);
      pji[0] = ptr[0];
      pji[2] = ptr[1];
      pji[1] = ptr[2];
      pji[3] = ptr[3];

      ptr += 4;

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  gets the derivative vector contribution from the suplied dense adjoints
///
/// \param[out]     jaa      size _naa_var.  The derivative vector for aa jastrow variables.
/// \param[in]      daa      size 4*_norb*_norb (simple = _norb*_norb).  The dense adjoint array.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::dense_aa_adj_to_der_vec(S * const jaa,
                                                           const S * const daa) {

  if (_simple) {

    for (int i = 0; i < _norb; i++)
    for (int j = i; j < _norb; j++) {
      S & val = jaa[ formic::cmpd_pair_index_incl(i, j, _norb) ];
      val += daa[i+_norb*j];
      val += daa[j+_norb*i];
    }

  } else {

    //std::fill(daa, daa + 4*_norb*_norb, formic::zero(S()));
    S * ptr = jaa;
    for (int j = 0; j < _norb; j++)
    for (int i = 0; i < j; i++) {

      const S * const pij = daa + 4*(i+_norb*j);
      ptr[0] += pij[0];
      ptr[1] += pij[1];
      ptr[2] += pij[2];
      ptr[3] += pij[3];

      //S * const pij = daa + 4*(i+_norb*j);
      //pij[0] = ptr[0];
      //pij[1] = ptr[1];
      //pij[2] = ptr[2];
      //pij[3] = ptr[3];

      const S * const pji = daa + 4*(j+_norb*i);
      ptr[0] += pji[0];
      ptr[1] += pji[2];
      ptr[2] += pji[1];
      ptr[3] += pji[3];

      //S * const pji = daa + 4*(j+_norb*i);
      //pji[0] = ptr[0];
      //pji[2] = ptr[1];
      //pji[1] = ptr[2];
      //pji[3] = ptr[3];

      ptr += 4;

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief builds the dense alpha-beta and beta-alpha jastrow factor arrays
///
/// \param[in]      jab      size _nab_var.  The compact array.
/// \param[out]     dab      size 4*_norb*_norb (simple = _norb*_norb).  On exit, the dense array in ab order.
/// \param[out]     dba      size 4*_norb*_norb (simple = _norb*_norb).  On exit, the dense array in ba order.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::fill_dense_ab(const S * const jab,
                                                       S * const dab,
                                                       S * const dba) {

  if (_simple) {

    if (_restricted) {

      for (int i = 0; i < _norb; i++)
      for (int j = i; j < _norb; j++) {
        const S val = jab[ formic::cmpd_pair_index_incl(i, j, _norb) ];
        dab[i+_norb*j] = val;
        dab[j+_norb*i] = val;
        dba[i+_norb*j] = val;
        dba[j+_norb*i] = val;
      }

    } else {

      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < _norb; i++) {
        const S val = jab[ i * _norb + j ];
        dab[i+_norb*j] = val;
        dba[j+_norb*i] = val;
      }

    }

  } else {

    for (int j = 0; j < _norb; j++)
    for (int i = 0; i < _norb; i++) {

      const S * const ptr = jab + 4*(i+_norb*j);

      S * const pij = dab + 4*(i+_norb*j);
      pij[0] = ptr[0];
      pij[1] = ptr[1];
      pij[2] = ptr[2];
      pij[3] = ptr[3];

      S * const pji = dba + 4*(j+_norb*i);
      pji[0] = ptr[0];
      pji[2] = ptr[1];
      pji[1] = ptr[2];
      pji[3] = ptr[3];

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief builds the dense alpha-beta and beta-alpha jastrow factor arrays
///
/// \param[out]     jab      size _nab_var.  The derivative vector for ab jastrow variables.
/// \param[in]      dab      size 4*_norb*_norb (simple = _norb*_norb).  The dense ab adjoint array
/// \param[in]      dba      size 4*_norb*_norb (simple = _norb*_norb).  The dense ba adjoint array
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::dense_ab_adj_to_der_vec(S * const jab,
                                                           const S * const dab,
                                                           const S * const dba) {

  if (_simple) {

    if (_restricted) {

      for (int i = 0; i < _norb; i++)
      for (int j = i; j < _norb; j++) {
        S & val = jab[ formic::cmpd_pair_index_incl(i, j, _norb) ];
        val += dab[i+_norb*j];
        val += dab[j+_norb*i];
        val += dba[i+_norb*j];
        val += dba[j+_norb*i];
        //const S val = jab[ formic::cmpd_pair_index_incl(i, j, _norb) ];
        //dab[i+_norb*j] = val;
        //dab[j+_norb*i] = val;
        //dba[i+_norb*j] = val;
        //dba[j+_norb*i] = val;
      }

    } else {

      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < _norb; i++) {
        S & val = jab[ i * _norb + j ];
        val += dab[i+_norb*j];
        val += dba[j+_norb*i];
        //const S val = jab[ i + _norb * j ];
        //dab[i+_norb*j] = val;
        //dba[j+_norb*i] = val;
      }

    }

  } else {

    for (int j = 0; j < _norb; j++)
    for (int i = 0; i < _norb; i++) {

      S * const ptr = jab + 4*(i+_norb*j);

      const S * const pij = dab + 4*(i+_norb*j);
      ptr[0] += pij[0];
      ptr[1] += pij[1];
      ptr[2] += pij[2];
      ptr[3] += pij[3];

      const S * const pji = dba + 4*(j+_norb*i);
      ptr[0] += pji[0];
      ptr[1] += pji[2];
      ptr[2] += pji[1];
      ptr[3] += pji[3];

      //const S * const ptr = jab + 4*(i+_norb*j);

      //S * const pij = dab + 4*(i+_norb*j);
      //pij[0] = ptr[0];
      //pij[1] = ptr[1];
      //pij[2] = ptr[2];
      //pij[3] = ptr[3];

      //S * const pji = dba + 4*(j+_norb*i);
      //pji[0] = ptr[0];
      //pji[2] = ptr[1];
      //pji[1] = ptr[2];
      //pji[3] = ptr[3];

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief constructs the coefficient object by reading names from an archive and finding the
///        jastrow factor arrays in the supplied wavefunction
///
/// \param[in]      userinp  user input options
/// \param[in,out]  wfn      the wavefunction holding the pairing matrix data
/// \param[in,out]  arch     an archive from which to read names
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
formic::fqmc::AllPairJast<S>::AllPairJast(const formic::InputBase & userinp, formic::Wavefunction<S> & wfn, formic::Archive & arch)
  : _norb(userinp.get<int>("nao")),
    _nap(userinp.get<int>("nalpha")),
    _nbp(userinp.get<int>("nbeta")),
    _nau(userinp.get<int>("nao") - userinp.get<int>("nalpha")),
    _nbu(userinp.get<int>("nao") - userinp.get<int>("nbeta")),
    _simple(userinp.get<bool>("simple_jastrow")),
    _restricted(userinp.get<bool>("restrict_jastrow")),
    _initialized(false),
    _config(userinp)
{

  if ( _restricted && !_simple )
    throw formic::Exception("the case of restricted, non-simple jastrow factors is not currently implemented");

  // get numbers of aa variables
  _naa_var = ( _simple ? ( _norb * ( _norb + 1 ) / 2 ) : ( 4 * _norb * ( _norb - 1 ) / 2 ) );

  // get number of ab variables
  if ( _restricted )
    _nab_var = _norb * ( _norb + 1 ) / 2; // restricted currently assumes that we are also simple
  else
    _nab_var = ( _simple ? 1 : 4 ) * _norb * _norb;

  // read in the name to give the coefficient object
  formic::mpi::read_and_bcast(arch, this->_name, "Failed to read name in AllPairJast constructor.");

  // read in names of the jastrow data
  std::string aa_name, bb_name, ab_name;
  formic::mpi::read_and_bcast(arch, aa_name, "Failed to read name of alpha-alpha jastrow factor data object in AllPairJast constructor.");
  if (!_restricted)
    formic::mpi::read_and_bcast(arch, bb_name, "Failed to read name of beta-beta jastrow factor data object in AllPairJast constructor.");
  formic::mpi::read_and_bcast(arch, ab_name, "Failed to read name of alpha-beta jastrow factor data object in AllPairJast constructor.");

  // get references to the vectors holding the jastrow data
  std::vector<S> & aa_vec = wfn.get_wf(aa_name)->data_vec();
  std::vector<S> & bb_vec = ( _restricted ? wfn.get_wf(aa_name)->data_vec() : wfn.get_wf(bb_name)->data_vec() );
  std::vector<S> & ab_vec = wfn.get_wf(ab_name)->data_vec();

  // ensure the jastrow factor arrays are the correct size
  if ( aa_vec.size() != _naa_var )
    throw formic::Exception("alpha-alpha jastrow factor array \"%s\" is the wrong size in constructor for AllPairJast object \"%s\"") % aa_name % this->_name;
  if ( !_restricted && bb_vec.size() != _naa_var )
    throw formic::Exception(" beta-beta  jastrow factor array \"%s\" is the wrong size in constructor for AllPairJast object \"%s\"") % bb_name % this->_name;
  if ( ab_vec.size() != _nab_var )
    throw formic::Exception("alpha-beta  jastrow factor array \"%s\" is the wrong size in constructor for AllPairJast object \"%s\"") % ab_name % this->_name;

  // set the jastrow array pointers
  _jfaa = &aa_vec.at(0);
  _jfbb = &bb_vec.at(0);
  _jfab = &ab_vec.at(0);

  // get whether the jastrow factor data is fixed
  _aa_fixed = wfn.get_wf(aa_name)->fixed();
  _bb_fixed = ( _restricted ? _aa_fixed : wfn.get_wf(bb_name)->fixed() );
  _ab_fixed = wfn.get_wf(ab_name)->fixed();

  // initialize dense jastrow factor arrays
  _djfaa.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  _djfbb.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  _djfab.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  _djfba.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  this->fill_dense_aa(_jfaa, &_djfaa.at(0));
  this->fill_dense_aa(_jfbb, &_djfbb.at(0));
  this->fill_dense_ab(_jfab, &_djfab.at(0), &_djfba.at(0));

  // initialize intermediate arrays

  _Jsi_a.assign(          _nap,  formic::zero(S()) );
  _Jsi_b.assign(          _nbp,  formic::zero(S()) );
  _Jsa_a.assign(          _nau,  formic::zero(S()) );
  _Jsa_b.assign(          _nbu,  formic::zero(S()) );
  _Jdij_aa.assign( _nap * _nap,  formic::zero(S()) );
  _Jdij_bb.assign( _nbp * _nbp,  formic::zero(S()) );
  _Jdij_ab.assign( _nap * _nbp,  formic::zero(S()) );
  _Jdab_aa.assign( _nau * _nau,  formic::zero(S()) );
  _Jdab_bb.assign( _nbu * _nbu,  formic::zero(S()) );
  _Jdab_ab.assign( _nau * _nbu,  formic::zero(S()) );
  _Jdia_aa.assign( _nap * _nau,  formic::zero(S()) );
  _Jdia_bb.assign( _nbp * _nbu,  formic::zero(S()) );
  _Jdia_ab.assign( _nap * _nbu,  formic::zero(S()) );
  _Jdai_aa.assign( _nau * _nap,  formic::zero(S()) );
  _Jdai_bb.assign( _nbu * _nbp,  formic::zero(S()) );
  _Jdai_ab.assign( _nau * _nbp,  formic::zero(S()) );

  _dzfaa.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  _dzfbb.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  _dzfab.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );
  _dzfba.assign( ( _simple ? 1 : 4 ) * _norb * _norb, formic::zero(S()) );

  _Zne = formic::zero(S());
  _Zsi_a.assign(          _nap,  formic::zero(S()) );
  _Zsi_b.assign(          _nbp,  formic::zero(S()) );
  _Zsa_a.assign(          _nau,  formic::zero(S()) );
  _Zsa_b.assign(          _nbu,  formic::zero(S()) );
  _Zdij_aa.assign( _nap * _nap,  formic::zero(S()) );
  _Zdij_bb.assign( _nbp * _nbp,  formic::zero(S()) );
  _Zdij_ab.assign( _nap * _nbp,  formic::zero(S()) );
  _Zdab_aa.assign( _nau * _nau,  formic::zero(S()) );
  _Zdab_bb.assign( _nbu * _nbu,  formic::zero(S()) );
  _Zdab_ab.assign( _nau * _nbu,  formic::zero(S()) );
  _Zdia_aa.assign( _nap * _nau,  formic::zero(S()) );
  _Zdia_bb.assign( _nbp * _nbu,  formic::zero(S()) );
  _Zdia_ab.assign( _nap * _nbu,  formic::zero(S()) );
  _Zdai_aa.assign( _nau * _nap,  formic::zero(S()) );
  _Zdai_bb.assign( _nbu * _nbp,  formic::zero(S()) );
  _Zdai_ab.assign( _nau * _nbp,  formic::zero(S()) );

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Initializes the object for the given lattice configuration.
///
/// \param[in]      config  the lattice configuration we want the object initialized for
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::initialize(const formic::Configuration & config) {

  // save the configuration
  _config = config;

  // initialize intermediates
  this->init_intermediates();

  // set the initialized flag to true
  _initialized = true;

}

template void formic::fqmc::AllPairJast<double>::initialize(const formic::Configuration &);
template void formic::fqmc::AllPairJast<std::complex<double> >::initialize(const formic::Configuration &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Updates the object so that the current configuration becomes the new configuration,
///        i.e. current = current + diff
///
/// \param[in]      diff    an object describing the change to the configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::update(const formic::ConfigDiff & diff) {


  // check that the object is initialized
  this->check_initialized("AllPairJast::update");

  // update the configuration
  _config.update(diff);

  // initialize intermediates
  this->init_intermediates();

}

template void formic::fqmc::AllPairJast<double>::update(const formic::ConfigDiff &);
template void formic::fqmc::AllPairJast<std::complex<double> >::update(const formic::ConfigDiff &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Returns the sum of the values in the jastrow factor exponent for the supplied
///        configuration, optionally including jastrow variables that are fixed.
///
/// \param[in]      config         the configuration to compute the sum for
/// \param[in]      jaa            the alpha-alpha jastrow factor variables to use
/// \param[in]      jbb            the beta-beta   jastrow factor variables to use
/// \param[in]      jab            the alpha-beta  jastrow factor variables to use
/// \param[in]      include_fixed  whether to include fixed jastrow variables
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
S formic::fqmc::AllPairJast<S>::get_exponent_sum(const formic::Configuration & config,
                                                 const S * const jaa,
                                                 const S * const jbb,
                                                 const S * const jab,
                                                 const bool include_fixed) {

  // initialize the sum of jastrow exponents
  S exp_sum = formic::zero(S());

  // simple jastrows
  if (_simple) {

    // get alpha and beta occupation vectors
    const int * const ao2s = &config.a_occ_to_site().at(0);
    const int * const bo2s = &config.b_occ_to_site().at(0);

    // alpha-alpha jastrows
    if ( include_fixed || !_aa_fixed ) {
      for (int i = 0; i < _nap; i++)
      for (int j = i; j < _nap; j++) {
        const int ii = std::min(ao2s[i], ao2s[j]);
        const int jj = std::max(ao2s[i], ao2s[j]);
        exp_sum += jaa[ formic::cmpd_pair_index_incl(ii, jj, _norb) ];
      }
    }

    // beta-beta jastrows
    if ( include_fixed || !_bb_fixed ) {
      for (int i = 0; i < _nbp; i++)
      for (int j = i; j < _nbp; j++) {
        const int ii = std::min(bo2s[i], bo2s[j]);
        const int jj = std::max(bo2s[i], bo2s[j]);
        exp_sum += jbb[ formic::cmpd_pair_index_incl(ii, jj, _norb) ];
      }
    }

    // alpha-beta jastrows
    if ( include_fixed || !_ab_fixed ) {
      if (_restricted) {
        for (int i = 0; i < _nap; i++)
        for (int j = 0; j < _nbp; j++) {
          const int ii = std::min(ao2s[i], bo2s[j]);
          const int jj = std::max(ao2s[i], bo2s[j]);
          exp_sum += jab[ formic::cmpd_pair_index_incl(ii, jj, _norb) ];
        }
      } else {
        for (int i = 0; i < _nap; i++)
        for (int j = 0; j < _nbp; j++) {
          const int ii = ao2s[i];
          const int jj = bo2s[j];
          exp_sum += jab[ ii * _norb + jj ];
        }
      }
    }

  // non-simple jastrows
  } else {

    // aa
    if ( include_fixed || !_aa_fixed ) {
      int ij = 0;
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < j; i++)
        exp_sum += jaa[ 4*(ij++) + config[2*i+0] + 2*config[2*j+0] ];
    }

    // bb
    if ( include_fixed || !_bb_fixed ) {
      int ij = 0;
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < j; i++)
        exp_sum += jbb[ 4*(ij++) + config[2*i+1] + 2*config[2*j+1] ];
    }

    // ab
    if ( include_fixed || !_ab_fixed ) {
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < _norb; i++)
        exp_sum += jab[ 4*(i+_norb*j) + config[2*i+0] + 2*config[2*j+1] ];
    }

  }

  // return the sum of exponents
  return exp_sum;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief returns the value of the coefficient for an arbitrary configuration
///
/// \param[in]      config  the lattice configuration for which we want the coefficient
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
formic::MantExp<S> formic::fqmc::AllPairJast<S>::get_coeff(const formic::Configuration & config) {

//  // initialize the sum of jastrow exponents
//  S exp_sum = formic::zero(S());
//
//  // coeff for simple jastrows
//  if (_simple) {
//
//    // get alpha and beta occupation vectors
//    const int * const ao2s = &config.a_occ_to_site().at(0);
//    const int * const bo2s = &config.b_occ_to_site().at(0);
//
//    // get convenient pointers to the jastrow data
//    const S * const jaa = _jfaa;
//    const S * const jbb = ( _restricted ? _jfaa : _jfbb );
//    const S * const jab = _jfab;
//
//    // alpha-alpha jastrows
//    for (int i = 0; i < _nap; i++)
//    for (int j = i; j < _nap; j++) {
//      const int ii = std::min(ao2s[i], ao2s[j]);
//      const int jj = std::max(ao2s[i], ao2s[j]);
//      exp_sum += jaa[ formic::cmpd_pair_index_incl(ii, jj, _norb) ];
//    }
//
//    // beta-beta jastrows
//    for (int i = 0; i < _nbp; i++)
//    for (int j = i; j < _nbp; j++) {
//      const int ii = std::min(bo2s[i], bo2s[j]);
//      const int jj = std::max(bo2s[i], bo2s[j]);
//      exp_sum += jbb[ formic::cmpd_pair_index_incl(ii, jj, _norb) ];
//    }
//
//    // alpha-beta jastrows
//    if (_restricted) {
//      for (int i = 0; i < _nap; i++)
//      for (int j = 0; j < _nbp; j++) {
//        const int ii = std::min(ao2s[i], bo2s[j]);
//        const int jj = std::max(ao2s[i], bo2s[j]);
//        exp_sum += jab[ formic::cmpd_pair_index_incl(ii, jj, _norb) ];
//      }
//    } else {
//      for (int i = 0; i < _nap; i++)
//      for (int j = 0; j < _nbp; j++) {
//        const int ii = ao2s[i];
//        const int jj = bo2s[j];
//        exp_sum += jab[ ii + _norb * jj ];
//      }
//    }
//
//  // coeff for non-simple jastrows
//  } else {
//
//    // aa
//    int ij = 0;
//    for (int j = 0; j < _norb; j++)
//    for (int i = 0; i < j; i++)
//      exp_sum += _jfaa[ 4*(ij++) + config[2*i+0] + 2*config[2*j+0] ];
//
//    // bb
//    ij = 0;
//    for (int j = 0; j < _norb; j++)
//    for (int i = 0; i < j; i++)
//      exp_sum += _jfbb[ 4*(ij++) + config[2*i+1] + 2*config[2*j+1] ];
//
//    // ab
//    for (int j = 0; j < _norb; j++)
//    for (int i = 0; i < _norb; i++)
//      exp_sum += _jfab[ 4*(i+_norb*j) + config[2*i+0] + 2*config[2*j+1] ];
//
//  }
//
//  // return the exponent of the sum of jastrow coefficients
//  return formic::MantExp<S>(std::exp(exp_sum));

  // return the exponent of the sum of jastrow exponents
  return formic::MantExp<S>( std::exp( this->get_exponent_sum(config, _jfaa, _jfbb, _jfab, true) ) );

}

template formic::MantExp<double> formic::fqmc::AllPairJast<double>::get_coeff(const formic::Configuration &);
template formic::MantExp<std::complex<double> > formic::fqmc::AllPairJast<std::complex<double> >::get_coeff(const formic::Configuration &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief returns the ratio of the coefficient at a new configuration (current + diff) to
///        that at the current configuration
///
/// \param[in]      diff     the differece between the new and current configurations
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
S formic::fqmc::AllPairJast<S>::get_ratio(const formic::ConfigDiff & diff) {

  // check that the object is initialized
  this->check_initialized("AllPairJast::get_ratio");

  // check that empty and full lists are the same size
  if ( diff.newly_full().size() != diff.newly_empty().size() )
    throw formic::Exception("different sized newly empty and newly full lists in AllPairJast::get_ratio");

  // check that the configuration change is possible
  formic::check_config_update(diff, _config, "AllPairJast::get_ratio");

  //formic::Configuration conf2 = _config;
  //conf2.update(diff);
  //return std::exp( this->get_exponent_sum(conf2, _jfaa, _jfbb, _jfab, true) - this->get_exponent_sum(_config, _jfaa, _jfbb, _jfab, true) );

  // no change
  if ( diff.newly_full().size() == 0 ) {

    return formic::unity(S());

  // single electron move
  } else if ( diff.newly_full().size() == 1 ) {

    const int I = diff.newly_empty().at(0);
    const int A = diff.newly_full().at(0);
    if ( I % 2 == 0 ) {
      const int i = _config.a_site_to_occ().at(I/2);
      const int a = _config.a_site_to_unocc().at(A/2);
      return _Jsi_a[i] * _Jsa_a[a] * _Jdia_aa[i+_nap*a];
    } else {
      const int i = _config.b_site_to_occ().at(I/2);
      const int a = _config.b_site_to_unocc().at(A/2);
      return _Jsi_b[i] * _Jsa_b[a] * _Jdia_bb[i+_nbp*a];
    }

  // double electron move
  } else if ( diff.newly_full().size() == 2 ) {

    int I = diff.newly_empty().at(0);
    int A = diff.newly_full().at(0);
    int J = diff.newly_empty().at(1);
    int B = diff.newly_full().at(1);

    if ( I % 2 != A % 2 || J % 2 != B % 2 )
      throw formic::Exception("spin mismatch in double move in AllPairJast::get_ratio");

    // aaaa
    if ( I % 2 == 0 && J % 2 == 0 ) {
      const int i = _config.a_site_to_occ().at(I/2);
      const int a = _config.a_site_to_unocc().at(A/2);
      const int j = _config.a_site_to_occ().at(J/2);
      const int b = _config.a_site_to_unocc().at(B/2);
      return   _Jsi_a[i] * _Jsi_a[j] * _Jsa_a[a] * _Jsa_a[b]
             * _Jdij_aa[i+_nap*j] * _Jdia_aa[i+_nap*a] * _Jdia_aa[i+_nap*b]
             * _Jdab_aa[a+_nau*b] * _Jdia_aa[j+_nap*a] * _Jdia_aa[j+_nap*b];
    // bbbb
    } else if ( I % 2 == 1 && J % 2 == 1 ) {
      const int i = _config.b_site_to_occ().at(I/2);
      const int a = _config.b_site_to_unocc().at(A/2);
      const int j = _config.b_site_to_occ().at(J/2);
      const int b = _config.b_site_to_unocc().at(B/2);
      return   _Jsi_b[i] * _Jsi_b[j] * _Jsa_b[a] * _Jsa_b[b]
             * _Jdij_bb[i+_nbp*j] * _Jdia_bb[i+_nbp*a] * _Jdia_bb[i+_nbp*b]
             * _Jdab_bb[a+_nbu*b] * _Jdia_bb[j+_nbp*a] * _Jdia_bb[j+_nbp*b];
    // aabb or bbaa
    } else {
      if ( I % 2 == 1 ) {
        assert( J % 2 == 0 );
        std::swap(I,J);
        std::swap(A,B);
      }
      const int i = _config.a_site_to_occ().at(I/2);
      const int a = _config.a_site_to_unocc().at(A/2);
      const int j = _config.b_site_to_occ().at(J/2);
      const int b = _config.b_site_to_unocc().at(B/2);
      return   _Jsi_a[i] * _Jsi_b[j] * _Jsa_a[a] * _Jsa_b[b]
             * _Jdij_ab[i+_nap*j] * _Jdia_aa[i+_nap*a] * _Jdia_ab[i+_nap*b]
             * _Jdab_ab[a+_nau*b] * _Jdai_ab[a+_nau*j] * _Jdia_bb[j+_nbp*b];
    }

  }

  throw formic::Exception("AllPairJast::get_ratio not implemented for 3 or more electron moves");

  return formic::unity(S());

}

template double formic::fqmc::AllPairJast<double>::get_ratio(const formic::ConfigDiff &);
template std::complex<double>  formic::fqmc::AllPairJast<std::complex<double> >::get_ratio(const formic::ConfigDiff &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief sets the pointers to the aa, bb, and ab variables inside an all-variable vector
///
/// \param[in]      vec          the all-variable vector, size == this->n_der()
/// \param[out]     aa_pos       a pointer to the aa jastrow position in the vector
/// \param[out]     bb_pos       a pointer to the bb jastrow position in the vector
/// \param[out]     ab_pos       a pointer to the ab jastrow position in the vector
/// \param[in]      ign_fixed    whether to ignore "fixed" variables
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::get_aa_bb_ab_position_pointers(S * const vec, S * & aa_pos, S * & bb_pos, S * & ab_pos, const bool ign_fixed) const {

  // aa variables are first
  aa_pos = vec;

  // bb variables are the same as aa variables if the jastrows are restricted, otherwise they come after
  bb_pos = aa_pos + ( _restricted || ( ign_fixed && _aa_fixed ) ? 0 : _naa_var );

  // last come the ab variables
  ab_pos = bb_pos + ( ign_fixed && _bb_fixed ? 0 : _naa_var );

}

template void formic::fqmc::AllPairJast<double>::get_aa_bb_ab_position_pointers(double * const vec, double * & aa_pos, double * & bb_pos, double * & ab_pos, const bool) const;
template void formic::fqmc::AllPairJast<std::complex<double> >::get_aa_bb_ab_position_pointers(std::complex<double>  * const vec, std::complex<double>  * & aa_pos, std::complex<double>  * & bb_pos, std::complex<double>  * & ab_pos, const bool) const;

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief sets the pointers to the aa, bb, and ab variables inside an all-variable vector
///
/// \param[in]      vec      the all-variable vector, size == this->n_der()
/// \param[out]     aa_pos   a pointer to the aa jastrow position in the vector
/// \param[out]     bb_pos   a pointer to the bb jastrow position in the vector
/// \param[out]     ab_pos   a pointer to the ab jastrow position in the vector
/// \param[in]      ign_fixed    whether to ignore "fixed" variables
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::get_aa_bb_ab_position_pointers(const S * const vec, const S * & aa_pos, const S * & bb_pos, const S * & ab_pos, const bool ign_fixed) const {

  // aa variables are first
  aa_pos = vec;

  // bb variables are the same as aa variables if the jastrows are restricted, otherwise they come after
  bb_pos = aa_pos + ( _restricted || ( ign_fixed && _aa_fixed ) ? 0 : _naa_var );

  // last come the ab variables
  ab_pos = bb_pos + ( ign_fixed && _bb_fixed ? 0 : _naa_var );

}

template void formic::fqmc::AllPairJast<double>::get_aa_bb_ab_position_pointers(const double * const vec, const double * & aa_pos, const double * & bb_pos, const double * & ab_pos, const bool) const;
template void formic::fqmc::AllPairJast<std::complex<double> >::get_aa_bb_ab_position_pointers(const std::complex<double>  * const vec, const std::complex<double>  * & aa_pos, const std::complex<double>  * & bb_pos, const std::complex<double>  * & ab_pos, const bool) const;

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief gives the derivative ratios for the current configuration
///
/// \param[out]     ratios     on exit, the derivative ratios. size == this->n_der()
///                            for ign_fixed == true, size == this->n_der()
///                            for ign_fixed == false, size == this->n_var()
/// \param[in]      ign_fixed  whether to ignore variables that are fixed
///
///   The ratio for variable i is     <n|Psi_i> / <n|Psi>
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::der_ratios(S * const ratios, const bool ign_fixed) {

  // check that the object is initialized
  this->check_initialized("AllPairJast::der_ratios");

  // get convenient pointers to the ratios
  S * aa_ratios;
  S * bb_ratios;
  S * ab_ratios;
  this->get_aa_bb_ab_position_pointers(ratios, aa_ratios, bb_ratios, ab_ratios, ign_fixed);

  // coeff for simple jastrows
  if (_simple) {

    // get alpha and beta occupation vectors
    const int * const ao2s = &_config.a_occ_to_site().at(0);
    const int * const bo2s = &_config.b_occ_to_site().at(0);

    // alpha-alpha jastrows
    if ( !_aa_fixed || !ign_fixed ) {
      std::fill(aa_ratios, aa_ratios + _naa_var, formic::zero(S()));
      for (int i = 0; i < _nap; i++)
      for (int j = i; j < _nap; j++) {
        const int ii = std::min(ao2s[i], ao2s[j]);
        const int jj = std::max(ao2s[i], ao2s[j]);
        aa_ratios[ formic::cmpd_pair_index_incl(ii, jj, _norb) ] += formic::unity(S());
      }
    }

    // beta-beta jastrows
    if ( !_bb_fixed || !ign_fixed ) {
      if (!_restricted)
        std::fill(bb_ratios, bb_ratios + _naa_var, formic::zero(S()));
      for (int i = 0; i < _nbp; i++)
      for (int j = i; j < _nbp; j++) {
        const int ii = std::min(bo2s[i], bo2s[j]);
        const int jj = std::max(bo2s[i], bo2s[j]);
        bb_ratios[ formic::cmpd_pair_index_incl(ii, jj, _norb) ] += formic::unity(S());
      }
    }

    // alpha-beta jastrows
    if ( !_ab_fixed || !ign_fixed ) {
      std::fill(ab_ratios, ab_ratios + _nab_var, formic::zero(S()));
      if (_restricted) {
        for (int i = 0; i < _nap; i++)
        for (int j = 0; j < _nbp; j++) {
          const int ii = std::min(ao2s[i], bo2s[j]);
          const int jj = std::max(ao2s[i], bo2s[j]);
          ab_ratios[ formic::cmpd_pair_index_incl(ii, jj, _norb) ] += formic::unity(S());
        }
      } else {
        for (int i = 0; i < _nap; i++)
        for (int j = 0; j < _nbp; j++) {
          const int ii = ao2s[i];
          const int jj = bo2s[j];
          ab_ratios[ ii * _norb + jj ] += formic::unity(S());
        }
      }
    }

  // coeff for non-simple jastrows
  } else {

    // aa
    if ( !_aa_fixed || !ign_fixed ) {
      std::fill(aa_ratios, aa_ratios + _naa_var, formic::zero(S()));
      int ij = 0;
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < j; i++)
        aa_ratios[ 4*(ij++) + _config[2*i+0] + 2*_config[2*j+0] ] = formic::unity(S());
    }

    // bb
    if ( !_bb_fixed || !ign_fixed ) {
      std::fill(bb_ratios, bb_ratios + _naa_var, formic::zero(S()));
      int ij = 0;
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < j; i++)
        bb_ratios[ 4*(ij++) + _config[2*i+1] + 2*_config[2*j+1] ] = formic::unity(S());
    }

    // ab
    if ( !_ab_fixed || !ign_fixed ) {
      std::fill(ab_ratios, ab_ratios + _nab_var, formic::zero(S()));
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < _norb; i++)
        ab_ratios[ 4*(i+_norb*j) + _config[2*i+0] + 2*_config[2*j+1] ] = formic::unity(S());
    }

  }

}

template void formic::fqmc::AllPairJast<double>::der_ratios(double * const, const bool ign_fixed);
template void formic::fqmc::AllPairJast<std::complex<double> >::der_ratios(std::complex<double>  * const, const bool ign_fixed);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief updates the coefficient object's variables
///
/// \param[in]      dfc      the derivative function coefficients.  size == this->n_der()
/// \param[in]      alpha    the constant to scale the update by
///
///  |Psi(x_1, x_2, ...)> + sum_{i=1} dfc_i |dPsi/dx_i>
///     --->  |Psi(x_1 + alpha * dfc_1, x_2 + alpha * dfc_2, ...)>
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::update_variables_detail(const S * const dfc, const S alpha) {

  // get convenient pointers to the data components
  const S * aa_update;
  const S * bb_update;
  const S * ab_update;
  this->get_aa_bb_ab_position_pointers(dfc, aa_update, bb_update, ab_update);

  // add the updates
  if ( !_aa_fixed                 ) formic::xaxpy(_naa_var, alpha, aa_update, 1, _jfaa, 1);
  if ( !_bb_fixed && !_restricted ) formic::xaxpy(_naa_var, alpha, bb_update, 1, _jfbb, 1);
  if ( !_ab_fixed                 ) formic::xaxpy(_nab_var, alpha, ab_update, 1, _jfab, 1);

//  // update for simple jastrows
//  if (_simple) {
//
//    // add the updates
//    if ( !_aa_fixed                 ) formic::xaxpy(_naa_var, alpha, aa_update, 1, _jfaa, 1);
//    if ( !_bb_fixed && !_restricted ) formic::xaxpy(_naa_var, alpha, bb_update, 1, _jfbb, 1);
//    if ( !_ab_fixed                 ) formic::xaxpy(_nab_var, alpha, ab_update, 1, _jfab, 1);
//
//  // update for non-simple jastrows
//  } else {
//
//    if (!_aa_fixed) formic::xaxpy(_naa_var, alpha, aa_update, 1, _jfaa, 1);
//    if (!_bb_fixed) formic::xaxpy(_naa_var, alpha, bb_update, 1, _jfbb, 1);
//    if (!_ab_fixed) formic::xaxpy(_nab_var, alpha, ab_update, 1, _jfab, 1);
//
//  }

  // fill the dense jastrow arrays with the new variables
  this->fill_dense_aa(_jfaa, &_djfaa.at(0));
  this->fill_dense_aa(_jfbb, &_djfbb.at(0));
  this->fill_dense_ab(_jfab, &_djfab.at(0), &_djfba.at(0));

  // clear the exact sampling vectors
  this->clear_exact_sample_vecs();

}

template void formic::fqmc::AllPairJast<double>::update_variables_detail(const double * const, const double);
template void formic::fqmc::AllPairJast<std::complex<double> >::update_variables_detail(const std::complex<double> * const, const std::complex<double>);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  combines the alpha-first and beta-first derivative vectors into the alpha-first
///
/// \param[in,out]  afirst   the derivatives with alpha first. size == this->n_der()
/// \param[in]      bfirst   the derivatives with  beta first. size == this->n_der()
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::add_beta_derivs_to_alpha_derivs(S * const afirst, const S * const bfirst) {

  // get convenient pointers to aa bb and ab vector positions
  S * af_aa_pos;
  S * af_bb_pos;
  S * af_ab_pos;
  this->get_aa_bb_ab_position_pointers(afirst, af_aa_pos, af_bb_pos, af_ab_pos);
  const S * bf_aa_pos;
  const S * bf_bb_pos;
  const S * bf_ab_pos;
  this->get_aa_bb_ab_position_pointers(bfirst, bf_bb_pos, bf_aa_pos, bf_ab_pos);

  // simple jastrows
  if (_simple) {

    if (_restricted) {
      formic::xaxpy(_naa_var, formic::unity(S()), bf_aa_pos, 1, af_aa_pos, 1);
      formic::xaxpy(_nab_var, formic::unity(S()), bf_ab_pos, 1, af_ab_pos, 1);
    } else {
      formic::xaxpy(_naa_var, formic::unity(S()), bf_aa_pos, 1, af_aa_pos, 1);
      formic::xaxpy(_naa_var, formic::unity(S()), bf_bb_pos, 1, af_bb_pos, 1);
      for (int j = 0; j < _norb; j++)
      for (int i = 0; i < _norb; i++)
        af_ab_pos[i*_norb+j] += bf_ab_pos[j*_norb+i];
    }

  // non-simple jastrows
  } else {

    throw formic::Exception("formic::fqmc::AllPairJast<S>::add_beta_derivs_to_alpha_derivs currently only implemented for simple jastrows");

    //// copy the data from the vector
    //int p = 0;
    //if (!_aa_fixed) formic::xcopy(_naa_var, vars + (p++)*_naa_var, 1, _jfaa, 1);
    //if (!_bb_fixed) formic::xcopy(_naa_var, vars + (p++)*_naa_var, 1, _jfbb, 1);
    //if (!_ab_fixed) formic::xcopy(_nab_var, vars + (p++)*_naa_var, 1, _jfab, 1);

  }

}

template void formic::fqmc::AllPairJast<double>::add_beta_derivs_to_alpha_derivs(double * const afirst, const double * const bfirst);
template void formic::fqmc::AllPairJast<std::complex<double> >::add_beta_derivs_to_alpha_derivs(std::complex<double>  * const afirst, const std::complex<double>  * const bfirst);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief copies the coefficient's variables to the supplied vector
///
/// \param[out]     vars     on exit, the variables.
/// \param[in]      all      whether to copy all variables or just non-fixed variables
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::to_vector_detail(S * const vars, const bool all) {

  // get convenient pointers to aa bb and ab vector positions
  S * aa_pos;
  S * bb_pos;
  S * ab_pos;
  this->get_aa_bb_ab_position_pointers(vars, aa_pos, bb_pos, ab_pos, !all);

  // copy the data to the vector
  if ( ( all || !_aa_fixed )                 ) formic::xcopy(_naa_var, _jfaa, 1, aa_pos, 1);
  if ( ( all || !_bb_fixed ) && !_restricted ) formic::xcopy(_naa_var, _jfbb, 1, bb_pos, 1);
  if ( ( all || !_ab_fixed )                 ) formic::xcopy(_nab_var, _jfab, 1, ab_pos, 1);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief copies the coefficient's variables to the supplied vector
///
/// \param[out]     vars     on exit, the variables. size == this->n_der()
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::to_vector(S * const vars) {

  this->to_vector_detail(vars, false);

}

template void formic::fqmc::AllPairJast<double>::to_vector(double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::to_vector(std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief copies the coefficient's variables (including "fixed" ones) to the supplied vector
///
/// \param[out]     vars     on exit, the variables. size == this->n_var()
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::all_to_vector(S * const vars) {

  this->to_vector_detail(vars, true);

}

template void formic::fqmc::AllPairJast<double>::all_to_vector(double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::all_to_vector(std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief copies the coefficient's variables from the supplied vector
///
/// \param[in]      vars     the variables.
/// \param[in]      all      whether to copy all or just the non-fixed variables
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::from_vector_detail(const S * const vars, const bool all) {

  // get convenient pointers to aa bb and ab vector positions
  const S * aa_pos;
  const S * bb_pos;
  const S * ab_pos;
  this->get_aa_bb_ab_position_pointers(vars, aa_pos, bb_pos, ab_pos, !all);

  // copy the data from the vector
  if ( ( all || !_aa_fixed )                 ) formic::xcopy(_naa_var, aa_pos, 1, _jfaa, 1);
  if ( ( all || !_bb_fixed ) && !_restricted ) formic::xcopy(_naa_var, bb_pos, 1, _jfbb, 1);
  if ( ( all || !_ab_fixed )                 ) formic::xcopy(_nab_var, ab_pos, 1, _jfab, 1);

  // fill the dense jastrow arrays with the new values
  this->fill_dense_aa(_jfaa, &_djfaa.at(0));
  this->fill_dense_aa(_jfbb, &_djfbb.at(0));
  this->fill_dense_ab(_jfab, &_djfab.at(0), &_djfba.at(0));

  // clear the exact sampling vectors
  this->clear_exact_sample_vecs();

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief copies the coefficient's variables from the supplied vector
///
/// \param[in]      vars     the variables. size == this->n_der()
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::from_vector(const S * const vars) {

  this->from_vector_detail(vars, false);

}

template void formic::fqmc::AllPairJast<double>::from_vector(const double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::from_vector(const std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief copies the coefficient's variables (including "fixed" ones) from the supplied vector
///
/// \param[in]      vars     the variables. size == this->n_var()
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::all_from_vector(const S * const vars) {

  this->from_vector_detail(vars, true);

}

template void formic::fqmc::AllPairJast<double>::all_from_vector(const double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::all_from_vector(const std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief prints the coefficient's variables to the output file
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::print_variables() {

  struct local_funcs {

    static void print_upper_triangle(const int n, const S * const lt_matrix) {
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < i; j++)
          formic::of << boost::format(" %14s") % ' ';
        for (int j = i; j < n; j++)
          formic::of << boost::format(" %s") % formic::format_number("%14.8f", lt_matrix[formic::cmpd_pair_index_incl(i, j, n)]);
        formic::of << std::endl;
      }
    }

    static void print_cstyle_matrix(const int n, const S * const matrix) {
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
          formic::of << boost::format(" %s") % formic::format_number("%14.8f", matrix[i*n+j]);
        formic::of << std::endl;
      }
    }

  };

  if ( formic::mpi::rank() == 0 ) {

    //formic::of << "printing alpha-alpha jastrow factor matrix:" << std::endl;
    //formic::of << boost::format("%3s ") % "";
    //for (int i = 0; i < _norb; i++)
    //  formic::of << boost::format("%5s%3i %12s") % "" % i % "";
    //formic::of << std::endl;
    //formic::of << std::endl;
    //for (int i = 0; i < _norb; i++) {
    //  for (int nj = 0; nj < 2; nj++)
    //  for (int ni = 0; ni < 2; ni++) {
    //    formic::of << boost::format("%3s ") % ( ni == 0 && nj == 0 ? (boost::format("%3i") % i).str() : std::string("") );
    //    for (int j = 0; j < _norb; j++) {
    //      if (j > i)
    //        formic::of << boost::format("%5s %i%i %s") % "" % ni % nj % formic::format_number("%12.6f", _djfaa[4*(i+_norb*j)+ni+2*nj]);
    //      else
    //        formic::of << boost::format("%5s %2s %12s") % "" % "" % "";
    //    }
    //    formic::of << std::endl;
    //  }
    //  formic::of << std::endl;
    //}

    // simple jastrows
    if (_simple) {

      // alpha-alpha or same-spin
      if (_restricted)
        formic::of << "printing same-spin jastrow factor matrix:" << std::endl;
      else
        formic::of << "printing alpha-alpha jastrow factor matrix:" << std::endl;
      local_funcs::print_upper_triangle(_norb, _jfaa);
      formic::of << std::endl;

      // beta-beta
      if (!_restricted) {
        formic::of << "printing beta-beta jastrow factor matrix:" << std::endl;
        local_funcs::print_upper_triangle(_norb, _jfbb);
        formic::of << std::endl;
      }

      // alpha-beta or opposite-spin
      if (_restricted) {
        formic::of << "printing opposite-spin jastrow factor matrix:" << std::endl;
        local_funcs::print_upper_triangle(_norb, _jfab);
      } else {
        formic::of << "printing alpha-beta jastrow factor matrix:" << std::endl;
        local_funcs::print_cstyle_matrix(_norb, _jfab);
      }
      formic::of << std::endl;

    // non-simple jastrows
    } else {

      formic::of << "printing alpha-alpha jastrow factor matrix:" << std::endl;
      for (int i = 0; i < _naa_var; i++)
        formic::of << boost::format("%s\n") % formic::format_number("%20.12f", _jfaa[i]);
      formic::of << std::endl;
      formic::of << "printing  beta-beta  jastrow factor matrix:" << std::endl;
      for (int i = 0; i < _naa_var; i++)
        formic::of << boost::format("%s\n") % formic::format_number("%20.12f", _jfbb[i]);
      formic::of << std::endl;
      formic::of << "printing alpha-beta  jastrow factor matrix:" << std::endl;
      for (int i = 0; i < _nab_var; i++)
        formic::of << boost::format("%s\n") % formic::format_number("%20.12f", _jfab[i]);
      formic::of << std::endl;

    }

  }

}

template void formic::fqmc::AllPairJast<double>::print_variables();
template void formic::fqmc::AllPairJast<std::complex<double> >::print_variables();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief initializes a single-index intermediate array
///
/// \param[out]     val_array     size n.  on exit, the single intermediate array
/// \param[in]      n             length of the intermediate array
/// \param[in]      do_exp        whether or not to exponentiate the array
/// \param[in]      ja            jastrow factor dense array when other spin is alpha
/// \param[in]      jb            jastrow factor dense array when other spin is beta
/// \param[in]      js            jastrow factor dense array for the same-site term
/// \param[in]      do_same_site  whether to explicitly add the same-site term
/// \param[in]      o2s           map from ordering in occ/unocc vector to lattice site
/// \param[in]      ao2s          map from ordering in alpha occ vector to lattice site
/// \param[in]      bo2s          map from ordering in beta  occ vector to lattice site
/// \param[in]      intmd_sign    sign to apply to the intermediate
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_single_intermed_simple(S * const val_array,
                                                               const int n,
                                                               const bool do_exp,
                                                               const S * const ja,
                                                               const S * const jb,
                                                               const S * const js,
                                                               const bool do_same_site,
                                                               const int * const o2s,
                                                               const int * const ao2s,
                                                               const int * const bo2s,
                                                               const double intmd_sign) {

  // prepare zero
  const S zz = formic::zero(S());

  // loop over the index of the changing site
  for (int i = 0; i < n; i++) {

    // get the dense array index offset for the current orbital
    const int offset = _norb * o2s[i];

    // get a temporary pointer for use in fast looping
    const S * jt;

    // initialize the value
    S x = zz;

    // contribution when other site is alpha
    jt = ja + offset;
    for (int j = 0; j < _nap; j++)
      x += jt[ ao2s[j] ];

    // contribution when other site is beta
    jt = jb + offset;
    for (int j = 0; j < _nbp; j++)
      x += jt[ bo2s[j] ];

    // add the same-site term if required
    if (do_same_site)
      x += js[ o2s[i] + offset ];

    // apply the sign
    x *= intmd_sign;

    // record total, exponentiating if requested
    val_array[i] = ( do_exp ? std::exp(x) : x );

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  expands the derivative adjoint of a single-index intermediate into the adjoints of the
///         dense jastrow factor arrays
///         It is assumed that the intermediate array was exponentiated.
///
/// \param[in]      adj_array     size n.  the single intermediate adjoint array
/// \param[in]      val_array     size n.  the single intermediate value array
/// \param[in]      n             length of the intermediate array
/// \param[in,out]  ja            jastrow factor adjoint dense array when other spin is alpha
/// \param[in,out]  jb            jastrow factor adjoint dense array when other spin is beta
/// \param[in,out]  js            jastrow factor adjoint dense array for the same-site term
/// \param[in]      do_same_site  whether the same-site term was used
/// \param[in]      o2s           map from ordering in occ/unocc vector to lattice site
/// \param[in]      ao2s          map from ordering in alpha occ vector to lattice site
/// \param[in]      bo2s          map from ordering in beta  occ vector to lattice site
/// \param[in]      intmd_sign    sign that was applied to the intermediate before exponentiating
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::expand_adj_single_intermed_simple(const S * const adj_array,
                                                                     const S * const val_array,
                                                                     const int n,
                                                                           S * const jaa,
                                                                           S * const jab,
                                                                           //S * const ja,
                                                                           //S * const jb,
                                                                           //S * const js,
                                                                     const bool do_same_site,
                                                                     const int * const o2s,
                                                                     const int * const ao2s,
                                                                     const int * const bo2s,
                                                                     const double intmd_sign) {

  // loop over the index of the changing site
  for (int i = 0; i < n; i++) {

    // get the adjoint times derivative value
    const S x = intmd_sign * val_array[i] * adj_array[i];

    // contribution when other site is alpha
    for (int j = 0; j < _nap; j++) {
      const int p = std::min(ao2s[j], o2s[i]);
      const int q = std::max(ao2s[j], o2s[i]);
      jaa[ formic::cmpd_pair_index_incl(p, q, _norb) ] += x;
    }

    // contribution when other site is beta
    if (_restricted) {
      for (int j = 0; j < _nap; j++) {
        const int p = std::min(bo2s[j], o2s[i]);
        const int q = std::max(bo2s[j], o2s[i]);
        jab[ formic::cmpd_pair_index_incl(p, q, _norb) ] += x;
      }
    } else {
      for (int j = 0; j < _nap; j++)
        jab[ _norb * o2s[i] + bo2s[j] ] += x;
    }

    // add the same-site term if required
    if (do_same_site)
      jaa[ formic::cmpd_pair_index_incl(o2s[i], o2s[i], _norb) ] += x;

  }

//  // loop over the index of the changing site
//  for (int i = 0; i < n; i++) {
//
//    // get the adjoint times derivative value
//    const S x = intmd_sign * val_array[i] * adj_array[i];
//
//    // get the dense array index offset for the current orbital
//    const int offset = _norb * o2s[i];
//
//    // get a temporary pointer for use in fast looping
//    S * jt;
//
//    // contribution when other site is alpha
//    jt = ja + offset;
//    for (int j = 0; j < _nap; j++)
//      jt[ ao2s[j] ] += x;
//
//    // contribution when other site is beta
//    jt = jb + offset;
//    for (int j = 0; j < _nbp; j++)
//      jt[ bo2s[j] ] += x;
//
//    // add the same-site term if required
//    if (do_same_site)
//      js[ o2s[i] + offset ] += x;
//
//  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief initializes a double-index intermediate array
///
/// \param[out]     val_array     size n1*n2.  on exit, the double intermediate array
/// \param[in]      n1            size of 1st dimension
/// \param[in]      n2            size of 2nd dimension
/// \param[in]      do_exp        whether or not to exponentiate the array
/// \param[in]      jd            jastrow factor dense array
/// \param[in]      o2s1          map from ordering in occ/unocc vector to lattice site
/// \param[in]      o2s2          map from ordering in occ/unocc vector to lattice site
/// \param[in]      intmd_sign    sign to apply to the intermediate
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_double_intermed_simple(S * const val_array,
                                                               const int n1,
                                                               const int n2,
                                                               const bool do_exp,
                                                               const S * const jd,
                                                               const int * const o2s1,
                                                               const int * const o2s2,
                                                               const double intmd_sign) {

  // loop over second index
  for (int j = 0; j < n2; j++) {

    // get temporary pointers for use in fast looping
    const S * const jt = jd + _norb * o2s2[j];
    S * const vt = val_array + n1 * j;

    // loop over the first index
    for (int i = 0; i < n1; i++) {

      // get jastrow element and apply the sign
      const S x = intmd_sign * jt[ o2s1[i] ];
      //const S x = intmd_sign * jd[ o2s1[i] + _norb * o2s2[j] ];

      // exponentiate if requested
      vt[i] = ( do_exp ? std::exp(x) : x );
      //val_array[i+n1*j] = ( do_exp ? std::exp(x) : x );

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  expands the derivative adjoint of a double-index intermediate into the adjoints of the
///         dense jastrow factor arrays
///         It is assumed that the intermediate array was exponentiated.
///
/// \param[in]      adj_array     size n.  the single intermediate adjoint array
/// \param[in]      val_array     size n.  the single intermediate value array
/// \param[in]      n1            size of 1st dimension
/// \param[in]      n2            size of 2nd dimension
/// \param[in,out]  jd            jastrow factor adjoint array
/// \param[in]      o2s1          map from ordering in occ/unocc vector to lattice site
/// \param[in]      o2s2          map from ordering in occ/unocc vector to lattice site
/// \param[in]      intmd_sign    sign to apply to the intermediate
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::expand_adj_double_intermed_simple_aa(const S * const adj_array,
                                                                        const S * const val_array,
                                                                        const int n1,
                                                                        const int n2,
                                                                              S * const jd,
                                                                        const int * const o2s1,
                                                                        const int * const o2s2,
                                                                        const double intmd_sign) {

  // loop over second index
  for (int j = 0; j < n2; j++) {

    // get temporary pointers for use in fast looping
    const S * const vt = val_array + n1 * j;
    const S * const at = adj_array + n1 * j;

    // loop over the first index
    for (int i = 0; i < n1; i++) {
      const int p = std::min(o2s2[j], o2s1[i]);
      const int q = std::max(o2s2[j], o2s1[i]);
      jd[ formic::cmpd_pair_index_incl(p, q, _norb) ] += intmd_sign * vt[i] * at[i];
    }

  }

//  // loop over second index
//  for (int j = 0; j < n2; j++) {
//
//    // get temporary pointers for use in fast looping
//    S * const jt = jd + _norb * o2s2[j];
//    const S * const vt = val_array + n1 * j;
//    const S * const at = adj_array + n1 * j;
//
//    // loop over the first index
//    for (int i = 0; i < n1; i++)
//      jt[ o2s1[i] ] += intmd_sign * vt[i] * at[i];
//
//  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  expands the derivative adjoint of a double-index intermediate into the adjoints of the
///         dense jastrow factor arrays
///         It is assumed that the intermediate array was exponentiated.
///
/// \param[in]      adj_array     size n.  the single intermediate adjoint array
/// \param[in]      val_array     size n.  the single intermediate value array
/// \param[in]      n1            size of 1st dimension
/// \param[in]      n2            size of 2nd dimension
/// \param[in,out]  jd            jastrow factor adjoint array
/// \param[in]      o2s1          map from ordering in occ/unocc vector to lattice site
/// \param[in]      o2s2          map from ordering in occ/unocc vector to lattice site
/// \param[in]      intmd_sign    sign to apply to the intermediate
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::expand_adj_double_intermed_simple_ab(const S * const adj_array,
                                                                        const S * const val_array,
                                                                        const int n1,
                                                                        const int n2,
                                                                              S * const jd,
                                                                        const int * const o2s1,
                                                                        const int * const o2s2,
                                                                        const double intmd_sign) {

  if (_restricted) {

    this->expand_adj_double_intermed_simple_aa(adj_array, val_array, n1, n2, jd, o2s1, o2s2, intmd_sign);

//    // loop over second index
//    for (int j = 0; j < n2; j++) {
//
//      // get temporary pointers for use in fast looping
//      const S * const vt = val_array + n1 * j;
//      const S * const at = adj_array + n1 * j;
//
//      // loop over the first index
//      for (int i = 0; i < n1; i++) {
//        const int p = std::min(o2s2[j], o2s1[i]);
//        const int q = std::max(o2s2[j], o2s1[i]);
//        jd[ formic::cmpd_pair_index_incl(p, q, _norb) ] += intmd_sign * vt[i] * at[i];
//      }
//
//    }

  } else {

    // loop over second index
    for (int j = 0; j < n2; j++) {

      // get temporary pointers for use in fast looping
      const S * const vt = val_array + n1 * j;
      const S * const at = adj_array + n1 * j;

      // loop over the first index
      for (int i = 0; i < n1; i++)
        jd[ _norb * o2s1[i] + o2s2[j] ] += intmd_sign * vt[i] * at[i];

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief initializes a single index intermediate array
///
/// \param[out]     val_array     size n.  on exit, the single intermediate array
/// \param[in]      n             length of the intermediate array
/// \param[in]      plus_offset   occupation offset for positive jastrow coefficient
/// \param[in]      minus_offset  occupation offset for negative jastrow coefficient
/// \param[in]      do_exp        whether or not to exponentiate the array
/// \param[in]      jaa           same-spin jastrow factor dense array
/// \param[in]      jba           opposite-spin jastrow factor dense array
/// \param[in]      o2s           map from ordering in occ/unocc vector to lattice site
/// \param[in]      begin_a       beginning position to the pointer to same-spin occupations
/// \param[in]      begin_b       beginning position to the pointer to opposite-spin occupations
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_single_intermed(S * const val_array,
                                                        const int n,
                                                        const int plus_offset,
                                                        const int minus_offset,
                                                        const bool do_exp,
                                                        const S * const jaa,
                                                        const S * const jba,
                                                        const int * const o2s,
                                                        const int * const begin_a,
                                                        const int * const begin_b) {

  // prepare some variables
  const S zz = formic::zero(S());
  const int * const end_a = begin_a + 2*_norb;
  const int * const end_b = begin_b + 2*_norb;
  int itmp, ptmp;
  const int * pocc;

  // loop over the index of the changing site
  for (int i = 0; i < n; i++) {

    // initialize the value
    S x = zz;

    // aa contribution, looping over aa jastrows touching our site
    itmp = 4*_norb*o2s[i];
    pocc = begin_a; // pointer to occupation of other site
    while (pocc != end_a) {
      ptmp = itmp + *pocc;
      x = x + jaa[ ptmp + plus_offset  ];
      x = x - jaa[ ptmp + minus_offset ];
      itmp += 4;
      pocc += 2;
    }

    // ab contribution, looping over ab jastrows touching our site
    itmp = 4*_norb*o2s[i];
    pocc = begin_b; // pointer to occupation of other site
    while (pocc != end_b) {
      ptmp = itmp + *pocc;
      x = x + jba[ ptmp + plus_offset ];
      x = x - jba[ ptmp + minus_offset ];
      itmp += 4;
      pocc += 2;
    }

    // record total, exponentiating if requested
    val_array[i] = ( do_exp ? std::exp(x) : x );

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief initializes a double index intermediate array
///
/// \param[out]     val_array     size ni * nj.  on exit, the double intermediate array
/// \param[in]      ni            range of the 1st index
/// \param[in]      nj            range of the 2nd index
/// \param[in]      ex_off_i      1st index occupation offset for excited determinant
/// \param[in]      rg_off_i      1st index occupation offset for regular determinant
/// \param[in]      ex_off_j      2nd index occupation offset for excited determinant
/// \param[in]      rg_off_j      2nd index occupation offset for regular determinant
/// \param[in]      do_exp        whether or not to exponentiate the array
/// \param[in]      jxx           jastrow factor dense array
/// \param[in]      ou2si         1st index occ/unocc to site mapping
/// \param[in]      ou2sj         2nd index occ/unocc to site mapping
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_double_intermed(S * const val_array,
                                                        const int ni,
                                                        const int nj,
                                                        const int ex_off_i,
                                                        const int rg_off_i,
                                                        const int ex_off_j,
                                                        const int rg_off_j,
                                                        const bool do_exp,
                                                        const S * const jxx,
                                                        const int * const ou2si,
                                                        const int * const ou2sj) {

  // prepare some variables
  int jtmp;
  const S * jf_ptr; // pointer to the coefficient array of a jastrow factor
  S x;

  // loop over the second index
  for (int j = 0; j < nj; j++) {

    jtmp = 4 * _norb * ou2sj[j];

    // loop over the first index
    for (int i = 0; i < ni; i++) {

      // get pointer to jastrow factor coefficient array
      jf_ptr = jxx + ( jtmp + 4 * ou2si[i] );

      // compute element
      x =     jf_ptr[ rg_off_i + rg_off_j ];
      x = x - jf_ptr[ ex_off_i + rg_off_j ];
      x = x - jf_ptr[ rg_off_i + ex_off_j ];
      x = x + jf_ptr[ ex_off_i + ex_off_j ];

      // record total, exponentiating if requested
      val_array[i+ni*j] = ( do_exp ? std::exp(x) : x );

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief initializes the supplied intermediate tensors
///
/// \param[in]      jaa      the dense alpha-alpha jastrow factor coefficient array
/// \param[in]      jbb      the dense  beta-beta  jastrow factor coefficient array
/// \param[in]      jab      the dense alpha-beta  jastrow factor coefficient array
/// \param[in]      jba      the dense  beta-alpha jastrow factor coefficient array
/// \param[out]     Jsi_a    intermediate
/// \param[out]     Jsi_b    intermediate
/// \param[out]     Jsa_a    intermediate
/// \param[out]     Jsa_b    intermediate
/// \param[out]     Jdij_aa  intermediate
/// \param[out]     Jdij_bb  intermediate
/// \param[out]     Jdij_ab  intermediate
/// \param[out]     Jdab_aa  intermediate
/// \param[out]     Jdab_bb  intermediate
/// \param[out]     Jdab_ab  intermediate
/// \param[out]     Jdia_aa  intermediate
/// \param[out]     Jdia_bb  intermediate
/// \param[out]     Jdia_ab  intermediate
/// \param[out]     Jdai_aa  intermediate
/// \param[out]     Jdai_bb  intermediate
/// \param[out]     Jdai_ab  intermediate
/// \param[in]      do_exp   whether to exponentiate the intermediate arrays
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_intermed_detail(const S * const jaa,
                                                        const S * const jbb,
                                                        const S * const jab,
                                                        const S * const jba,
                                                              S * const Jsi_a,
                                                              S * const Jsi_b,
                                                              S * const Jsa_a,
                                                              S * const Jsa_b,
                                                              S * const Jdij_aa,
                                                              S * const Jdij_bb,
                                                              S * const Jdij_ab,
                                                              S * const Jdab_aa,
                                                              S * const Jdab_bb,
                                                              S * const Jdab_ab,
                                                              S * const Jdia_aa,
                                                              S * const Jdia_bb,
                                                              S * const Jdia_ab,
                                                              S * const Jdai_aa,
                                                              S * const Jdai_bb,
                                                              S * const Jdai_ab,
                                                        const bool do_exp) {

  // get pointers to the occ/unocc to site mappings
  const int * const ao2s = &_config.a_occ_to_site().at(0);
  const int * const au2s = &_config.a_unocc_to_site().at(0);
  const int * const bo2s = &_config.b_occ_to_site().at(0);
  const int * const bu2s = &_config.b_unocc_to_site().at(0);

  // get pointer to the beginning of the configuration occupation number vector
  const int * const cfg_ptr = &_config.fock_config()[0];

  // simple jastrows
  if (_simple) {

    // prepare alpha single index intermediates
    this->init_single_intermed_simple(&Jsi_a[0], _nap, do_exp, jaa, jba, jaa, false, ao2s, ao2s, bo2s, -1.0);
    this->init_single_intermed_simple(&Jsa_a[0], _nau, do_exp, jaa, jba, jaa,  true, au2s, ao2s, bo2s,  1.0);

    // prepare beta single index intermediates
    this->init_single_intermed_simple(&Jsi_b[0], _nbp, do_exp, jab, jbb, jbb, false, bo2s, ao2s, bo2s, -1.0);
    this->init_single_intermed_simple(&Jsa_b[0], _nbu, do_exp, jab, jbb, jbb,  true, bu2s, ao2s, bo2s,  1.0);

    // prepare alpha-alpha double index intermediates
    this->init_double_intermed_simple(&Jdij_aa[0], _nap, _nap, do_exp, jaa, ao2s, ao2s,  1.0);
    this->init_double_intermed_simple(&Jdia_aa[0], _nap, _nau, do_exp, jaa, ao2s, au2s, -1.0);
    this->init_double_intermed_simple(&Jdai_aa[0], _nau, _nap, do_exp, jaa, au2s, ao2s, -1.0);
    this->init_double_intermed_simple(&Jdab_aa[0], _nau, _nau, do_exp, jaa, au2s, au2s,  1.0);

    // prepare beta-beta double index intermediates
    this->init_double_intermed_simple(&Jdij_bb[0], _nbp, _nbp, do_exp, jbb, bo2s, bo2s,  1.0);
    this->init_double_intermed_simple(&Jdia_bb[0], _nbp, _nbu, do_exp, jbb, bo2s, bu2s, -1.0);
    this->init_double_intermed_simple(&Jdai_bb[0], _nbu, _nbp, do_exp, jbb, bu2s, bo2s, -1.0);
    this->init_double_intermed_simple(&Jdab_bb[0], _nbu, _nbu, do_exp, jbb, bu2s, bu2s,  1.0);

    // prepare alpha-beta double index intermediates
    this->init_double_intermed_simple(&Jdij_ab[0], _nap, _nbp, do_exp, jab, ao2s, bo2s,  1.0);
    this->init_double_intermed_simple(&Jdia_ab[0], _nap, _nbu, do_exp, jab, ao2s, bu2s, -1.0);
    this->init_double_intermed_simple(&Jdai_ab[0], _nau, _nbp, do_exp, jab, au2s, bo2s, -1.0);
    this->init_double_intermed_simple(&Jdab_ab[0], _nau, _nbu, do_exp, jab, au2s, bu2s,  1.0);

  // non-simple jastrows
  } else {

    // prepare alpha single index intermediates
    this->init_single_intermed(&Jsi_a[0], _nap, 0, 2, do_exp, jaa, jba, ao2s, cfg_ptr+0, cfg_ptr+1);
    this->init_single_intermed(&Jsa_a[0], _nau, 2, 0, do_exp, jaa, jba, au2s, cfg_ptr+0, cfg_ptr+1);

    // prepare beta single index intermediates
    this->init_single_intermed(&Jsi_b[0], _nbp, 0, 2, do_exp, jbb, jab, bo2s, cfg_ptr+1, cfg_ptr+0);
    this->init_single_intermed(&Jsa_b[0], _nbu, 2, 0, do_exp, jbb, jab, bu2s, cfg_ptr+1, cfg_ptr+0);

    // prepare alpha-alpha double index intermediates
    this->init_double_intermed(&Jdij_aa[0], _nap, _nap, 0, 1, 0, 2, do_exp, jaa, ao2s, ao2s);
    this->init_double_intermed(&Jdia_aa[0], _nap, _nau, 0, 1, 2, 0, do_exp, jaa, ao2s, au2s);
    this->init_double_intermed(&Jdai_aa[0], _nau, _nap, 1, 0, 0, 2, do_exp, jaa, au2s, ao2s);
    this->init_double_intermed(&Jdab_aa[0], _nau, _nau, 1, 0, 2, 0, do_exp, jaa, au2s, au2s);

    // prepare beta-beta double index intermediates
    this->init_double_intermed(&Jdij_bb[0], _nbp, _nbp, 0, 1, 0, 2, do_exp, jbb, bo2s, bo2s);
    this->init_double_intermed(&Jdia_bb[0], _nbp, _nbu, 0, 1, 2, 0, do_exp, jbb, bo2s, bu2s);
    this->init_double_intermed(&Jdai_bb[0], _nbu, _nbp, 1, 0, 0, 2, do_exp, jbb, bu2s, bo2s);
    this->init_double_intermed(&Jdab_bb[0], _nbu, _nbu, 1, 0, 2, 0, do_exp, jbb, bu2s, bu2s);

    // prepare alpha-beta double index intermediates
    this->init_double_intermed(&Jdij_ab[0], _nap, _nbp, 0, 1, 0, 2, do_exp, jab, ao2s, bo2s);
    this->init_double_intermed(&Jdia_ab[0], _nap, _nbu, 0, 1, 2, 0, do_exp, jab, ao2s, bu2s);
    this->init_double_intermed(&Jdai_ab[0], _nau, _nbp, 1, 0, 0, 2, do_exp, jab, au2s, bo2s);
    this->init_double_intermed(&Jdab_ab[0], _nau, _nbu, 1, 0, 2, 0, do_exp, jab, au2s, bu2s);

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief initializes the intermediate tensors
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_intermediates() {

  this->init_intermed_detail(&_djfaa.at(0),
                             &_djfbb.at(0),
                             &_djfab.at(0),
                             &_djfba.at(0),
                             &_Jsi_a.at(0),
                             &_Jsi_b.at(0),
                             &_Jsa_a.at(0),
                             &_Jsa_b.at(0),
                             &_Jdij_aa.at(0),
                             &_Jdij_bb.at(0),
                             &_Jdij_ab.at(0),
                             &_Jdab_aa.at(0),
                             &_Jdab_bb.at(0),
                             &_Jdab_ab.at(0),
                             &_Jdia_aa.at(0),
                             &_Jdia_bb.at(0),
                             &_Jdia_ab.at(0),
                             &_Jdai_aa.at(0),
                             &_Jdai_bb.at(0),
                             &_Jdai_ab.at(0),
                             true);

}

template void formic::fqmc::AllPairJast<double>::init_intermediates();
template void formic::fqmc::AllPairJast<std::complex<double> >::init_intermediates();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  computes the derivative contributions from the supplied intermediate adjoints
///         for the aaaa case
///
/// \param[in]      ao2s        the alpha occ-to-site index array
/// \param[in]      bo2s        the alpha occ-to-site index array
/// \param[in]      au2s        the alpha unocc-to-site index array
/// \param[in]      bu2s        the alpha unocc-to-site index array
/// \param[in]      Jsi_a       intermediate
/// \param[in]      Jsa_a       intermediate
/// \param[in]      Jdij_aa     intermediate
/// \param[in]      Jdab_aa     intermediate
/// \param[in]      Jdia_aa     intermediate
/// \param[in]      Jdai_aa     intermediate
/// \param[in]      Asi_a       adjoint of intermediate
/// \param[in]      Asa_a       adjoint of intermediate
/// \param[in]      Adij_aa     adjoint of intermediate
/// \param[in]      Adab_aa     adjoint of intermediate
/// \param[in]      Adia_aa     adjoint of intermediate
/// \param[in]      Adai_aa     adjoint of intermediate
/// \param[in,out]  der_vec_aa  vector of derivatives of local energy for aa jastrow variables, size _naa_var
/// \param[in,out]  der_vec_ab  vector of derivatives of local energy for ab jastrow variables, size _nab_var
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::intermed_adj_to_der_vec_aa(const int * const ao2s,
                                                              const int * const bo2s,
                                                              const int * const au2s,
                                                              const int * const bu2s,
                                                              const S * const Jsi_a,
                                                              const S * const Jsa_a,
                                                              const S * const Jdij_aa,
                                                              const S * const Jdab_aa,
                                                              const S * const Jdia_aa,
                                                              const S * const Jdai_aa,
                                                              const S * const Asi_a,
                                                              const S * const Asa_a,
                                                              const S * const Adij_aa,
                                                              const S * const Adab_aa,
                                                              const S * const Adia_aa,
                                                              const S * const Adai_aa,
                                                                    S * const der_vec_aa,
                                                                    S * const der_vec_ab) {

  //// zero the dense jastrow arrays
  //std::fill(_dzfaa.begin(), _dzfaa.end(), formic::zero(S()));
  //std::fill(_dzfba.begin(), _dzfba.end(), formic::zero(S()));
  //std::fill(_dzfab.begin(), _dzfab.end(), formic::zero(S()));

  // get pointer to the beginning of the configuration occupation number vector
  const int * const cfg_ptr = &_config.fock_config()[0];

  // simple jastrows
  if (_simple) {

    // prepare alpha single index intermediates
    this->expand_adj_single_intermed_simple(Asi_a, Jsi_a, _nap, der_vec_aa, der_vec_ab, false, ao2s, ao2s, bo2s, -1.0);
    this->expand_adj_single_intermed_simple(Asa_a, Jsa_a, _nau, der_vec_aa, der_vec_ab,  true, au2s, ao2s, bo2s,  1.0);

    //this->expand_adj_single_intermed_simple(Asi_a, Jsi_a, _nap, &_dzfaa[0], &_dzfba[0], &_dzfaa[0], false, ao2s, ao2s, bo2s, -1.0);
    //this->expand_adj_single_intermed_simple(Asa_a, Jsa_a, _nau, &_dzfaa[0], &_dzfba[0], &_dzfaa[0],  true, au2s, ao2s, bo2s,  1.0);

    // prepare alpha-alpha double index intermediates
    this->expand_adj_double_intermed_simple_aa(Adij_aa, Jdij_aa, _nap, _nap, der_vec_aa, ao2s, ao2s,  1.0);
    this->expand_adj_double_intermed_simple_aa(Adia_aa, Jdia_aa, _nap, _nau, der_vec_aa, ao2s, au2s, -1.0);
    this->expand_adj_double_intermed_simple_aa(Adai_aa, Jdai_aa, _nau, _nap, der_vec_aa, au2s, ao2s, -1.0);
    this->expand_adj_double_intermed_simple_aa(Adab_aa, Jdab_aa, _nau, _nau, der_vec_aa, au2s, au2s,  1.0);

    //this->expand_adj_double_intermed_simple(Adij_aa, Jdij_aa, _nap, _nap, &_dzfaa[0], ao2s, ao2s,  1.0);
    //this->expand_adj_double_intermed_simple(Adia_aa, Jdia_aa, _nap, _nau, &_dzfaa[0], ao2s, au2s, -1.0);
    //this->expand_adj_double_intermed_simple(Adai_aa, Jdai_aa, _nau, _nap, &_dzfaa[0], au2s, ao2s, -1.0);
    //this->expand_adj_double_intermed_simple(Adab_aa, Jdab_aa, _nau, _nau, &_dzfaa[0], au2s, au2s,  1.0);

  // non-simple jastrows
  } else {

    throw formic::Exception("formic::fqmc::AllPairJast<S>::intermed_adj_to_der_vec_aa currently only implemented for simple jastrows");

//    // prepare alpha single index intermediates
//    this->init_single_intermed(&Jsi_a[0], _nap, 0, 2, do_exp, jaa, jba, ao2s, cfg_ptr+0, cfg_ptr+1);
//    this->init_single_intermed(&Jsa_a[0], _nau, 2, 0, do_exp, jaa, jba, au2s, cfg_ptr+0, cfg_ptr+1);
//
//    // prepare beta single index intermediates
//    this->init_single_intermed(&Jsi_b[0], _nbp, 0, 2, do_exp, jbb, jab, bo2s, cfg_ptr+1, cfg_ptr+0);
//    this->init_single_intermed(&Jsa_b[0], _nbu, 2, 0, do_exp, jbb, jab, bu2s, cfg_ptr+1, cfg_ptr+0);
//
//    // prepare alpha-alpha double index intermediates
//    this->init_double_intermed(&Jdij_aa[0], _nap, _nap, 0, 1, 0, 2, do_exp, jaa, ao2s, ao2s);
//    this->init_double_intermed(&Jdia_aa[0], _nap, _nau, 0, 1, 2, 0, do_exp, jaa, ao2s, au2s);
//    this->init_double_intermed(&Jdai_aa[0], _nau, _nap, 1, 0, 0, 2, do_exp, jaa, au2s, ao2s);
//    this->init_double_intermed(&Jdab_aa[0], _nau, _nau, 1, 0, 2, 0, do_exp, jaa, au2s, au2s);
//
//    // prepare beta-beta double index intermediates
//    this->init_double_intermed(&Jdij_bb[0], _nbp, _nbp, 0, 1, 0, 2, do_exp, jbb, bo2s, bo2s);
//    this->init_double_intermed(&Jdia_bb[0], _nbp, _nbu, 0, 1, 2, 0, do_exp, jbb, bo2s, bu2s);
//    this->init_double_intermed(&Jdai_bb[0], _nbu, _nbp, 1, 0, 0, 2, do_exp, jbb, bu2s, bo2s);
//    this->init_double_intermed(&Jdab_bb[0], _nbu, _nbu, 1, 0, 2, 0, do_exp, jbb, bu2s, bu2s);
//
//    // prepare alpha-beta double index intermediates
//    this->init_double_intermed(&Jdij_ab[0], _nap, _nbp, 0, 1, 0, 2, do_exp, jab, ao2s, bo2s);
//    this->init_double_intermed(&Jdia_ab[0], _nap, _nbu, 0, 1, 2, 0, do_exp, jab, ao2s, bu2s);
//    this->init_double_intermed(&Jdai_ab[0], _nau, _nbp, 1, 0, 0, 2, do_exp, jab, au2s, bo2s);
//    this->init_double_intermed(&Jdab_ab[0], _nau, _nbu, 1, 0, 2, 0, do_exp, jab, au2s, bu2s);

  }

  //// add contributions from dense adjoints to the derivative vectors
  //this->dense_aa_adj_to_der_vec(der_vec_aa, &_dzfaa[0]);
  //this->dense_ab_adj_to_der_vec(der_vec_ab, &_dzfab[0], &_dzfba[0]);

}

template void formic::fqmc::AllPairJast<double>::intermed_adj_to_der_vec_aa(const int * const ao2s,
                                                                    const int * const bo2s,
                                                                    const int * const au2s,
                                                                    const int * const bu2s,
                                                                    const double * const Jsi_a,
                                                                    const double * const Jsa_a,
                                                                    const double * const Jdij_aa,
                                                                    const double * const Jdab_aa,
                                                                    const double * const Jdia_aa,
                                                                    const double * const Jdai_aa,
                                                                    const double * const Asi_a,
                                                                    const double * const Asa_a,
                                                                    const double * const Adij_aa,
                                                                    const double * const Adab_aa,
                                                                    const double * const Adia_aa,
                                                                    const double * const Adai_aa,
                                                                          double * const der_vec_aa,
                                                                          double * const der_vec_ab);

template void formic::fqmc::AllPairJast<std::complex<double> >::intermed_adj_to_der_vec_aa(const int * const ao2s,
                                                                    const int * const bo2s,
                                                                    const int * const au2s,
                                                                    const int * const bu2s,
                                                                    const std::complex<double> * const Jsi_a,
                                                                    const std::complex<double> * const Jsa_a,
                                                                    const std::complex<double> * const Jdij_aa,
                                                                    const std::complex<double> * const Jdab_aa,
                                                                    const std::complex<double> * const Jdia_aa,
                                                                    const std::complex<double> * const Jdai_aa,
                                                                    const std::complex<double> * const Asi_a,
                                                                    const std::complex<double> * const Asa_a,
                                                                    const std::complex<double> * const Adij_aa,
                                                                    const std::complex<double> * const Adab_aa,
                                                                    const std::complex<double> * const Adia_aa,
                                                                    const std::complex<double> * const Adai_aa,
                                                                          std::complex<double> * const der_vec_aa,
                                                                          std::complex<double> * const der_vec_ab);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  computes the derivative contributions from the supplied intermediate adjoints
///         for the aaaa case
///
/// \param[in]      ao2s        the alpha occ-to-site index array
/// \param[in]      bo2s        the alpha occ-to-site index array
/// \param[in]      au2s        the alpha unocc-to-site index array
/// \param[in]      bu2s        the alpha unocc-to-site index array
/// \param[in]      Jsi_a       intermediate
/// \param[in]      Jsa_a       intermediate
/// \param[in]      Jdij_aa     intermediate
/// \param[in]      Jdab_aa     intermediate
/// \param[in]      Jdia_aa     intermediate
/// \param[in]      Jdai_aa     intermediate
/// \param[in]      Asi_a       adjoint of intermediate
/// \param[in]      Asa_a       adjoint of intermediate
/// \param[in]      Adij_aa     adjoint of intermediate
/// \param[in]      Adab_aa     adjoint of intermediate
/// \param[in]      Adia_aa     adjoint of intermediate
/// \param[in]      Adai_aa     adjoint of intermediate
/// \param[in,out]  der_vec_aa  vector of derivatives of local energy for aa jastrow variables, size _naa_var
/// \param[in,out]  der_vec_ab  vector of derivatives of local energy for ab jastrow variables, size _nab_var
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::intermed_adj_to_der_vec_ab(const int * const ao2s,
                                                              const int * const bo2s,
                                                              const int * const au2s,
                                                              const int * const bu2s,
                                                              const S * const Jsi_a,
                                                              const S * const Jsa_a,
                                                              const S * const Jsj_b,
                                                              const S * const Jsb_b,
                                                              const S * const Jdia_aa,
                                                              const S * const Jdij_ab,
                                                              const S * const Jdib_ab,
                                                              const S * const Jdaj_ab,
                                                              const S * const Jdab_ab,
                                                              const S * const Jdjb_bb,
                                                              const S * const Asi_a,
                                                              const S * const Asa_a,
                                                              const S * const Asj_b,
                                                              const S * const Asb_b,
                                                              const S * const Adia_aa,
                                                              const S * const Adij_ab,
                                                              const S * const Adib_ab,
                                                              const S * const Adaj_ab,
                                                              const S * const Adab_ab,
                                                              const S * const Adjb_bb,
                                                                    S * const der_vec_aa_a,
                                                                    S * const der_vec_ab_a,
                                                                    S * const der_vec_aa_b,
                                                                    S * const der_vec_ab_b) {

  // simple jastrows
  if (_simple) {

    // alpha single index
    this->expand_adj_single_intermed_simple(Asi_a, Jsi_a, _nap, der_vec_aa_a, der_vec_ab_a, false, ao2s, ao2s, bo2s, -1.0);
    this->expand_adj_single_intermed_simple(Asa_a, Jsa_a, _nau, der_vec_aa_a, der_vec_ab_a,  true, au2s, ao2s, bo2s,  1.0);

    // beta single index
    this->expand_adj_single_intermed_simple(Asj_b, Jsj_b, _nap, der_vec_aa_b, der_vec_ab_b, false, bo2s, bo2s, ao2s, -1.0);
    this->expand_adj_single_intermed_simple(Asb_b, Jsb_b, _nau, der_vec_aa_b, der_vec_ab_b,  true, bu2s, bo2s, ao2s,  1.0);

    // alpha-alpha double index
    this->expand_adj_double_intermed_simple_aa(Adia_aa, Jdia_aa, _nap, _nau, der_vec_aa_a, ao2s, au2s, -1.0);

    // beta-beta double index
    this->expand_adj_double_intermed_simple_aa(Adjb_bb, Jdjb_bb, _nap, _nau, der_vec_aa_b, bo2s, bu2s, -1.0);

    // alpha-beta double index
    this->expand_adj_double_intermed_simple_ab(Adij_ab, Jdij_ab, _nap, _nap, der_vec_ab_a, ao2s, bo2s,  1.0);
    this->expand_adj_double_intermed_simple_ab(Adib_ab, Jdib_ab, _nap, _nau, der_vec_ab_a, ao2s, bu2s, -1.0);
    this->expand_adj_double_intermed_simple_ab(Adaj_ab, Jdaj_ab, _nau, _nap, der_vec_ab_a, au2s, bo2s, -1.0);
    this->expand_adj_double_intermed_simple_ab(Adab_ab, Jdab_ab, _nau, _nau, der_vec_ab_a, au2s, bu2s,  1.0);

  // non-simple jastrows
  } else {

    throw formic::Exception("formic::fqmc::AllPairJast<S>::intermed_adj_to_der_vec_ab currently only implemented for simple jastrows");

  }

}

template void formic::fqmc::AllPairJast<double>::intermed_adj_to_der_vec_ab(const int * const ao2s,
                                                                       const int * const bo2s,
                                                                       const int * const au2s,
                                                                       const int * const bu2s,
                                                                       const double * const Jsi_a,
                                                                       const double * const Jsa_a,
                                                                       const double * const Jsj_b,
                                                                       const double * const Jsb_b,
                                                                       const double * const Jdia_aa,
                                                                       const double * const Jdij_ab,
                                                                       const double * const Jdib_ab,
                                                                       const double * const Jdaj_ab,
                                                                       const double * const Jdab_ab,
                                                                       const double * const Jdjb_bb,
                                                                       const double * const Asi_a,
                                                                       const double * const Asa_a,
                                                                       const double * const Asj_b,
                                                                       const double * const Asb_b,
                                                                       const double * const Adia_aa,
                                                                       const double * const Adij_ab,
                                                                       const double * const Adib_ab,
                                                                       const double * const Adaj_ab,
                                                                       const double * const Adab_ab,
                                                                       const double * const Adjb_bb,
                                                                             double * const der_vec_aa_a,
                                                                             double * const der_vec_ab_a,
                                                                             double * const der_vec_aa_b,
                                                                             double * const der_vec_ab_b);

template void formic::fqmc::AllPairJast<std::complex<double> >::intermed_adj_to_der_vec_ab(const int * const ao2s,
                                                                       const int * const bo2s,
                                                                       const int * const au2s,
                                                                       const int * const bu2s,
                                                                       const std::complex<double>  * const Jsi_a,
                                                                       const std::complex<double>  * const Jsa_a,
                                                                       const std::complex<double>  * const Jsj_b,
                                                                       const std::complex<double>  * const Jsb_b,
                                                                       const std::complex<double>  * const Jdia_aa,
                                                                       const std::complex<double>  * const Jdij_ab,
                                                                       const std::complex<double>  * const Jdib_ab,
                                                                       const std::complex<double>  * const Jdaj_ab,
                                                                       const std::complex<double>  * const Jdab_ab,
                                                                       const std::complex<double>  * const Jdjb_bb,
                                                                       const std::complex<double>  * const Asi_a,
                                                                       const std::complex<double>  * const Asa_a,
                                                                       const std::complex<double>  * const Asj_b,
                                                                       const std::complex<double>  * const Asb_b,
                                                                       const std::complex<double>  * const Adia_aa,
                                                                       const std::complex<double>  * const Adij_ab,
                                                                       const std::complex<double>  * const Adib_ab,
                                                                       const std::complex<double>  * const Adaj_ab,
                                                                       const std::complex<double>  * const Adab_ab,
                                                                       const std::complex<double>  * const Adjb_bb,
                                                                             std::complex<double>  * const der_vec_aa_a,
                                                                             std::complex<double>  * const der_vec_ab_a,
                                                                             std::complex<double>  * const der_vec_aa_b,
                                                                             std::complex<double>  * const der_vec_ab_b);

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_transform_oei(int ni,
                                                         int na,
                                                         double *   Jsi,
                                                         double *   Jsa,
                                                         double *  Jdia,
                                                         double *   Tia)
{
  FORMIC_FC_GLOBAL_(apj_fortran_transform_oei,APJ_FORTRAN_TRANSFORM_OEI)
    (&ni, &na, Jsi, Jsa, Jdia, Tia);
}

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_transform_oei(int ni,
                                                         int na,
                                                         std::complex<double> *   Jsi,
                                                         std::complex<double> *   Jsa,
                                                         std::complex<double> *  Jdia,
                                                         std::complex<double> *   Tia)
{
  throw formic::Exception("no implementation of AllPairJast::jastrow_transform_oei for complex numbers");
}

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_transform_tei(int ni,
                                                         int na,
                                                         int nj,
                                                         int nb,
                                                         double *   Jsi,
                                                         double *   Jsa,
                                                         double *   Jsj,
                                                         double *   Jsb,
                                                         double *  Jdia,
                                                         double *  Jdij,
                                                         double *  Jdib,
                                                         double *  Jdaj,
                                                         double *  Jdab,
                                                         double *  Jdjb,
                                                         double * Tiajb)
{
  FORMIC_FC_GLOBAL_(apj_fortran_transform_tei,APJ_FORTRAN_TRANSFORM_TEI)
    (&ni, &na, &nj, &nb, Jsi, Jsa, Jsj, Jsb, Jdia, Jdij, Jdib, Jdaj, Jdab, Jdjb, Tiajb);
}

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_transform_tei(int ni,
                                                         int na,
                                                         int nj,
                                                         int nb,
                                                         std::complex<double> *   Jsi,
                                                         std::complex<double> *   Jsa,
                                                         std::complex<double> *   Jsj,
                                                         std::complex<double> *   Jsb,
                                                         std::complex<double> *  Jdia,
                                                         std::complex<double> *  Jdij,
                                                         std::complex<double> *  Jdib,
                                                         std::complex<double> *  Jdaj,
                                                         std::complex<double> *  Jdab,
                                                         std::complex<double> *  Jdjb,
                                                         std::complex<double> * Tiajb)
{
  throw formic::Exception("no implementation of AllPairJast::jastrow_transform_tei for complex numbers");
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief 
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::transform_ints_aa_aaaa(S * const Tia, S * const Tiajb) {

  // one electron integrals
  this->jastrow_transform_oei(_nap,
                              _nau,
                              &_Jsi_a.at(0), 
                              &_Jsa_a.at(0), 
                              &_Jdia_aa.at(0), 
                              Tia);

  // two electron integrals
  this->jastrow_transform_tei(_nap,
                              _nau,
                              _nap,
                              _nau,
                              &_Jsi_a.at(0), 
                              &_Jsa_a.at(0), 
                              &_Jsi_a.at(0), 
                              &_Jsa_a.at(0), 
                              &_Jdia_aa.at(0), 
                              &_Jdij_aa.at(0), 
                              &_Jdia_aa.at(0), 
                              &_Jdai_aa.at(0), 
                              &_Jdab_aa.at(0), 
                              &_Jdia_aa.at(0), 
                              Tiajb);

}

template void formic::fqmc::AllPairJast<double>::transform_ints_aa_aaaa(double * const, double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::transform_ints_aa_aaaa(std::complex<double>  * const, std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief 
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::transform_ints_bb_bbbb(S * const Tia, S * const Tiajb) {

  // one electron integrals
  this->jastrow_transform_oei(_nbp,
                              _nbu,
                              &_Jsi_b.at(0), 
                              &_Jsa_b.at(0), 
                              &_Jdia_bb.at(0), 
                              Tia);

  // two electron integrals
  this->jastrow_transform_tei(_nbp,
                              _nbu,
                              _nbp,
                              _nbu,
                              &_Jsi_b.at(0), 
                              &_Jsa_b.at(0), 
                              &_Jsi_b.at(0), 
                              &_Jsa_b.at(0), 
                              &_Jdia_bb.at(0), 
                              &_Jdij_bb.at(0), 
                              &_Jdia_bb.at(0), 
                              &_Jdai_bb.at(0), 
                              &_Jdab_bb.at(0), 
                              &_Jdia_bb.at(0), 
                              Tiajb);

}

template void formic::fqmc::AllPairJast<double>::transform_ints_bb_bbbb(double * const, double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::transform_ints_bb_bbbb(std::complex<double>  * const, std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief 
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::transform_ints_aabb(S * const Tiajb) {

  // two electron integrals
  this->jastrow_transform_tei(_nap,
                              _nau,
                              _nbp,
                              _nbu,
                              &_Jsi_a.at(0), 
                              &_Jsa_a.at(0), 
                              &_Jsi_b.at(0), 
                              &_Jsa_b.at(0), 
                              &_Jdia_aa.at(0), 
                              &_Jdij_ab.at(0), 
                              &_Jdia_ab.at(0), 
                              &_Jdai_ab.at(0), 
                              &_Jdab_ab.at(0), 
                              &_Jdia_bb.at(0), 
                              Tiajb);

}

template void formic::fqmc::AllPairJast<double>::transform_ints_aabb(double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::transform_ints_aabb(std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Returns the identity contraction with a derivative ratio vector.  Note that the first
///        element of the vector corresponds to the no-derivative term.
///
/// \param[in]      dfc      the derivative function coefficients.  size == 1+this->n_der()
///
/// \return  sum_i <n|I|Psi_i> dfc_i / <n|Psi>
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
S formic::fqmc::AllPairJast<S>::der_ic(const S * const dfc) {

  // non-derivative term
  S retval = dfc[0];

  // get pointers to different jastrow derivative coefficients
  const S * dfc_aa;
  const S * dfc_bb;
  const S * dfc_ab;
  this->get_aa_bb_ab_position_pointers(dfc+1, dfc_aa, dfc_bb, dfc_ab);
//  if (_simple) {
//    this->get_aa_bb_ab_position_pointers(dfc+1, dfc_aa, dfc_bb, dfc_ab);
//  } else {
//    dfc_aa = dfc+1;
//    dfc_bb = dfc_aa + ( _aa_fixed ? 0 : _naa_var );
//    dfc_ab = dfc_bb + ( _bb_fixed ? 0 : _naa_var );
//  }

  // add the sum of non-fixed jastrow exponents
  retval += this->get_exponent_sum(_config, dfc_aa, dfc_bb, dfc_ab, false);

//  int p = 0;
//
//  // aa
//  if (!_aa_fixed) {
//    const S * const dfc_aa = dfc + 1 + (p++)*_naa_var;
//    int ij = 0;
//    for (int j = 0; j < _norb; j++)
//    for (int i = 0; i < j; i++)
//      retval += dfc_aa[ 4*(ij++) + _config[2*i+0] + 2*_config[2*j+0] ];
//  }
//
//  // bb
//  if (!_bb_fixed) {
//    const S * const dfc_bb = dfc + 1 + (p++)*_naa_var;
//    int ij = 0;
//    for (int j = 0; j < _norb; j++)
//    for (int i = 0; i < j; i++)
//      retval += dfc_bb[ 4*(ij++) + _config[2*i+1] + 2*_config[2*j+1] ];
//  }
//
//  // ab
//  if (!_ab_fixed) {
//    const S * const dfc_ab = dfc + 1 + (p++)*_naa_var;
//    for (int j = 0; j < _norb; j++)
//    for (int i = 0; i < _norb; i++)
//      retval += dfc_ab[ 4*(i+_norb*j) + _config[2*i+0] + 2*_config[2*j+1] ];
//  }

  // return the overall contraction
  return retval;

}

template double formic::fqmc::AllPairJast<double>::der_ic(const double * const);
template std::complex<double>  formic::fqmc::AllPairJast<std::complex<double> >::der_ic(const std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Prepares the intermediates needed to perform the integral transformation needed
///        for Hamiltonian derivative contraction.
///
/// \param[in]  dfc      the derivative function coefficients.  size == 1+this->n_der()
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::init_der_intermediates(const S * const dfc) {

  // get pointers to the different jastrows' derivative coefficients
  const S * dfc_aa;
  const S * dfc_bb;
  const S * dfc_ab;
  this->get_aa_bb_ab_position_pointers(dfc+1, dfc_aa, dfc_bb, dfc_ab);

  // prepare dense jastrow arrays from the derivative coefficients
  {

    // aa
    if (!_aa_fixed)
      this->fill_dense_aa(dfc_aa, &_dzfaa.at(0));
    else
      std::fill(_dzfaa.begin(), _dzfaa.end(), formic::zero(S()));

    // bb
    if (!_bb_fixed)
      this->fill_dense_aa(dfc_bb, &_dzfbb.at(0));
    else
      std::fill(_dzfbb.begin(), _dzfbb.end(), formic::zero(S()));

    // ab and ba
    if (!_ab_fixed) {
      this->fill_dense_ab(dfc_ab, &_dzfab.at(0), &_dzfba.at(0));
    } else {
      std::fill(_dzfab.begin(), _dzfab.end(), formic::zero(S()));
      std::fill(_dzfba.begin(), _dzfba.end(), formic::zero(S()));
    }

  }

//  // construct the dense arrays from the derivative function coefficients
//  {
//
//    int p = 0;
//
//    // aa
//    if (!_aa_fixed)
//      this->fill_dense_aa(dfc + 1 + (p++)*_naa_var, &_dzfaa.at(0));
//    else
//      std::fill(_dzfaa.begin(), _dzfaa.end(), formic::zero(S()));
//
//    // bb
//    if (!_bb_fixed)
//      this->fill_dense_aa(dfc + 1 + (p++)*_naa_var, &_dzfbb.at(0));
//    else
//      std::fill(_dzfbb.begin(), _dzfbb.end(), formic::zero(S()));
//
//    // ab and ba
//    if (!_ab_fixed) {
//      this->fill_dense_ab(dfc + 1 + (p++)*_naa_var, &_dzfab.at(0), &_dzfba.at(0));
//    } else {
//      std::fill(_dzfab.begin(), _dzfab.end(), formic::zero(S()));
//      std::fill(_dzfba.begin(), _dzfba.end(), formic::zero(S()));
//    }
//
//  }

  // initialize the derivative coefficient intermediate arrays
  this->init_intermed_detail(&_dzfaa.at(0),
                             &_dzfbb.at(0),
                             &_dzfab.at(0),
                             &_dzfba.at(0),
                             &_Zsi_a.at(0),
                             &_Zsi_b.at(0),
                             &_Zsa_a.at(0),
                             &_Zsa_b.at(0),
                             &_Zdij_aa.at(0),
                             &_Zdij_bb.at(0),
                             &_Zdij_ab.at(0),
                             &_Zdab_aa.at(0),
                             &_Zdab_bb.at(0),
                             &_Zdab_ab.at(0),
                             &_Zdia_aa.at(0),
                             &_Zdia_bb.at(0),
                             &_Zdia_ab.at(0),
                             &_Zdai_aa.at(0),
                             &_Zdai_bb.at(0),
                             &_Zdai_ab.at(0),
                             false);

  // initialize the non-excitation contraction
  _Zne = this->der_ic(dfc) - dfc[0];

}

template void formic::fqmc::AllPairJast<double>::init_der_intermediates(const double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::init_der_intermediates(const std::complex<double>  * const);
















template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_der_trans_oei(int ni,
                                                         int na,
                                                         double     Jne,
                                                         double *   Jsi,
                                                         double *   Jsa,
                                                         double *  Jdia,
                                                         double *   Tia)
{
  FORMIC_FC_GLOBAL_(apj_fortran_der_trans_oei,APJ_FORTRAN_DER_TRANS_OEI)
    (&ni, &na, &Jne, Jsi, Jsa, Jdia, Tia);
}

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_der_trans_oei(int ni,
                                                         int na,
                                                         std::complex<double>     Jne,
                                                         std::complex<double> *   Jsi,
                                                         std::complex<double> *   Jsa,
                                                         std::complex<double> *  Jdia,
                                                         std::complex<double> *   Tia)
{
  throw formic::Exception("no implementation of AllPairJast::jastrow_der_trans_oei for complex numbers");
}

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_der_trans_tei(int ni,
                                                         int na,
                                                         int nj,
                                                         int nb,
                                                         double     Jne,
                                                         double *   Jsi,
                                                         double *   Jsa,
                                                         double *   Jsj,
                                                         double *   Jsb,
                                                         double *  Jdia,
                                                         double *  Jdij,
                                                         double *  Jdib,
                                                         double *  Jdaj,
                                                         double *  Jdab,
                                                         double *  Jdjb,
                                                         double * Tiajb)
{
  FORMIC_FC_GLOBAL_(apj_fortran_der_trans_tei,APJ_FORTRAN_DER_TRANS_TEI)
    (&ni, &na, &nj, &nb, &Jne, Jsi, Jsa, Jsj, Jsb, Jdia, Jdij, Jdib, Jdaj, Jdab, Jdjb, Tiajb);
}

template<class S>
void formic::fqmc::AllPairJast<S>::jastrow_der_trans_tei(int ni,
                                                         int na,
                                                         int nj,
                                                         int nb,
                                                         std::complex<double>     Jne,
                                                         std::complex<double> *   Jsi,
                                                         std::complex<double> *   Jsa,
                                                         std::complex<double> *   Jsj,
                                                         std::complex<double> *   Jsb,
                                                         std::complex<double> *  Jdia,
                                                         std::complex<double> *  Jdij,
                                                         std::complex<double> *  Jdib,
                                                         std::complex<double> *  Jdaj,
                                                         std::complex<double> *  Jdab,
                                                         std::complex<double> *  Jdjb,
                                                         std::complex<double> * Tiajb)
{
  throw formic::Exception("no implementation of AllPairJast::jastrow_der_trans_tei for complex numbers");
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief 
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::der_trans_ints_aa_aaaa(S * const Tia, S * const Tiajb) {

  // one electron integrals
  this->jastrow_der_trans_oei(_nap,
                              _nau,
                              _Zne,
                              &_Zsi_a.at(0), 
                              &_Zsa_a.at(0), 
                              &_Zdia_aa.at(0), 
                              Tia);

  // two electron integrals
  this->jastrow_der_trans_tei(_nap,
                              _nau,
                              _nap,
                              _nau,
                              _Zne,
                              &_Zsi_a.at(0), 
                              &_Zsa_a.at(0), 
                              &_Zsi_a.at(0), 
                              &_Zsa_a.at(0), 
                              &_Zdia_aa.at(0), 
                              &_Zdij_aa.at(0), 
                              &_Zdia_aa.at(0), 
                              &_Zdai_aa.at(0), 
                              &_Zdab_aa.at(0), 
                              &_Zdia_aa.at(0), 
                              Tiajb);

}

template void formic::fqmc::AllPairJast<double>::der_trans_ints_aa_aaaa(double * const, double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::der_trans_ints_aa_aaaa(std::complex<double>  * const, std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief 
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::der_trans_ints_bb_bbbb(S * const Tia, S * const Tiajb) {

  // one electron integrals
  this->jastrow_der_trans_oei(_nbp,
                              _nbu,
                              _Zne,
                              &_Zsi_b.at(0), 
                              &_Zsa_b.at(0), 
                              &_Zdia_bb.at(0), 
                              Tia);

  // two electron integrals
  this->jastrow_der_trans_tei(_nbp,
                              _nbu,
                              _nbp,
                              _nbu,
                              _Zne,
                              &_Zsi_b.at(0), 
                              &_Zsa_b.at(0), 
                              &_Zsi_b.at(0), 
                              &_Zsa_b.at(0), 
                              &_Zdia_bb.at(0), 
                              &_Zdij_bb.at(0), 
                              &_Zdia_bb.at(0), 
                              &_Zdai_bb.at(0), 
                              &_Zdab_bb.at(0), 
                              &_Zdia_bb.at(0), 
                              Tiajb);

}

template void formic::fqmc::AllPairJast<double>::der_trans_ints_bb_bbbb(double * const, double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::der_trans_ints_bb_bbbb(std::complex<double>  * const, std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief 
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template<class S>
void formic::fqmc::AllPairJast<S>::der_trans_ints_aabb(S * const Tiajb) {

  // two electron integrals
  this->jastrow_der_trans_tei(_nap,
                              _nau,
                              _nbp,
                              _nbu,
                              _Zne,
                              &_Zsi_a.at(0), 
                              &_Zsa_a.at(0), 
                              &_Zsi_b.at(0), 
                              &_Zsa_b.at(0), 
                              &_Zdia_aa.at(0), 
                              &_Zdij_ab.at(0), 
                              &_Zdia_ab.at(0), 
                              &_Zdai_ab.at(0), 
                              &_Zdab_ab.at(0), 
                              &_Zdia_bb.at(0), 
                              Tiajb);

}

template void formic::fqmc::AllPairJast<double>::der_trans_ints_aabb(double * const);
template void formic::fqmc::AllPairJast<std::complex<double> >::der_trans_ints_aabb(std::complex<double>  * const);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  copy constructor is private and intended to be used by the clone function
///
/// \param[in]      other    the object to copy from
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S>
formic::fqmc::AllPairJast<S>::AllPairJast(const formic::fqmc::AllPairJast<S> & other)
  : formic::CoeffObj<S>(other),
    _norb(other._norb),
    _nap(other._nap),
    _nbp(other._nbp),
    _nau(other._nau),
    _nbu(other._nbu),
    _naa_var(other._naa_var),
    _nab_var(other._nab_var),
    _jfaa(other._jfaa),
    _jfbb(other._jfbb),
    _jfab(other._jfab),
    _djfaa(other._djfaa),
    _djfbb(other._djfbb),
    _djfab(other._djfab),
    _djfba(other._djfba),
    _simple(other._simple),
    _restricted(other._restricted),
    _aa_fixed(other._aa_fixed),
    _bb_fixed(other._bb_fixed),
    _ab_fixed(other._ab_fixed),
    _initialized(other._initialized),
    _config(other._config),
    _jfaa_data(other._jfaa_data),
    _jfbb_data(other._jfbb_data),
    _jfab_data(other._jfab_data),
    _Jsi_a(other._Jsi_a),
    _Jsi_b(other._Jsi_b),
    _Jsa_a(other._Jsa_a),
    _Jsa_b(other._Jsa_b),
    _Jdij_aa(other._Jdij_aa),
    _Jdij_bb(other._Jdij_bb),
    _Jdij_ab(other._Jdij_ab),
    _Jdab_aa(other._Jdab_aa),
    _Jdab_bb(other._Jdab_bb),
    _Jdab_ab(other._Jdab_ab),
    _Jdia_aa(other._Jdia_aa),
    _Jdia_bb(other._Jdia_bb),
    _Jdia_ab(other._Jdia_ab),
    _Jdai_aa(other._Jdai_aa),
    _Jdai_bb(other._Jdai_bb),
    _Jdai_ab(other._Jdai_ab),
    _dzfaa(other._dzfaa),
    _dzfbb(other._dzfbb),
    _dzfab(other._dzfab),
    _dzfba(other._dzfba),
    _Zne(other._Zne),
    _Zsi_a(other._Zsi_a),
    _Zsi_b(other._Zsi_b),
    _Zsa_a(other._Zsa_a),
    _Zsa_b(other._Zsa_b),
    _Zdij_aa(other._Zdij_aa),
    _Zdij_bb(other._Zdij_bb),
    _Zdij_ab(other._Zdij_ab),
    _Zdab_aa(other._Zdab_aa),
    _Zdab_bb(other._Zdab_bb),
    _Zdab_ab(other._Zdab_ab),
    _Zdia_aa(other._Zdia_aa),
    _Zdia_bb(other._Zdia_bb),
    _Zdia_ab(other._Zdia_ab),
    _Zdai_aa(other._Zdai_aa),
    _Zdai_bb(other._Zdai_bb),
    _Zdai_ab(other._Zdai_ab)
{

  // prepare the jf arrays (necessary if other's arrays point to the wavefunction)
  _jfaa_data.resize(_naa_var);
  _jfbb_data.resize(_naa_var);
  _jfab_data.resize(_nab_var);
  std::copy(other._jfaa, other._jfaa + _naa_var, _jfaa_data.begin());
  std::copy(other._jfbb, other._jfbb + _naa_var, _jfbb_data.begin());
  std::copy(other._jfab, other._jfab + _nab_var, _jfab_data.begin());
  _jfaa = &_jfaa_data.at(0);
  _jfbb = &_jfbb_data.at(0);
  _jfab = &_jfab_data.at(0);

}

template formic::fqmc::AllPairJast<double>::AllPairJast(const formic::fqmc::AllPairJast<double> &);
template formic::fqmc::AllPairJast<std::complex<double> >::AllPairJast(const formic::fqmc::AllPairJast<std::complex<double> > &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief clones the object
///
/// \return  a shared pointer to a new clone of the object
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S>
boost::shared_ptr<formic::CoeffObj<S> > formic::fqmc::AllPairJast<S>::clone() const {

  boost::shared_ptr<formic::CoeffObj<S> > retval( new formic::fqmc::AllPairJast<S>(*this) );
  return retval;

}

template boost::shared_ptr<formic::CoeffObj<double> > formic::fqmc::AllPairJast<double>::clone() const;
template boost::shared_ptr<formic::CoeffObj<std::complex<double> > > formic::fqmc::AllPairJast<std::complex<double> >::clone() const;
























///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "AllPairJast" command, which adds an all pairs jastrow
///          coefficient object to the wavefunction
///
///          The "AllPairJast" command takes the following arguments:
///            1. The first argument is the name to give the object.
///            2. The second argument is the name of the aa jastrow data array.
///            3. The third argument is the name of the bb jastrow data array.
///               Omit this argument if the jastrows are set to be restricted.
///            4. The final argument is the name of the ab jastrow data array.
///
/// \param[in,out]   arch     archive to read command arguments from
/// \param[in,out]   data     data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_AllPairJast_command_function(formic::Archive & arch, formic::AllData & data) {

  typedef std::complex<double> scdoub;

  // read in and add the restricted geminal power coefficient object to the wavefunction
  if (data.userinp.get<std::string>("scalar_type") == "real") {
    boost::shared_ptr<formic::CoeffObj<double> > c( new formic::fqmc::AllPairJast<double>(data.userinp, data.wd, arch) );
    data.wd.add_coeff_obj(c);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    boost::shared_ptr<formic::CoeffObj<scdoub> > c( new formic::fqmc::AllPairJast<scdoub>(data.userinp, data.wc, arch) );
    data.wc.add_coeff_obj(c);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" when processing AllPairJast command.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Expands the jastrow factors from a smaller orbital set and adds them to the jastrow
///          factors of a larger orbital set.
///
/// \param[in]       userinp  the user's input options
/// \param[in,out]   wfn      the wavefunction object holding the jastrow factor data
/// \param[in,out]   arch     archive to read command arguments from
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S>
static void formic_ExpandAPJ_command_function_detail(const formic::InputBase & userinp,
                                                     formic::Wavefunction<S> & wfn,
                                                     formic::Archive & arch) {

  // variables to describe the two sets of jastrow factors
  int lg_norb, sm_norb; // number of orbitals
  std::string lg_aa_name, lg_bb_name, lg_ab_name, sm_aa_name, sm_bb_name, sm_ab_name;

  // read in the number of orbitals and the jastrow factor names for the small set
  formic::mpi::read_and_bcast(arch, sm_norb, "Failed to read sm_norb in formic_ExpandAPJ_command_function_detail.");
  formic::mpi::read_and_bcast(arch, sm_aa_name, "Failed to read sm_aa_name in formic_ExpandAPJ_command_function_detail.");
  formic::mpi::read_and_bcast(arch, sm_bb_name, "Failed to read sm_bb_name in formic_ExpandAPJ_command_function_detail.");
  formic::mpi::read_and_bcast(arch, sm_ab_name, "Failed to read sm_ab_name in formic_ExpandAPJ_command_function_detail.");

  // read in the number of orbitals and the jastrow factor names for the large set
  formic::mpi::read_and_bcast(arch, lg_norb, "Failed to read lg_norb in formic_ExpandAPJ_command_function_detail.");
  formic::mpi::read_and_bcast(arch, lg_aa_name, "Failed to read lg_aa_name in formic_ExpandAPJ_command_function_detail.");
  formic::mpi::read_and_bcast(arch, lg_bb_name, "Failed to read lg_bb_name in formic_ExpandAPJ_command_function_detail.");
  formic::mpi::read_and_bcast(arch, lg_ab_name, "Failed to read lg_ab_name in formic_ExpandAPJ_command_function_detail.");

  // read in the index mapping
  std::vector<int> imap;
  formic::mpi::read_and_bcast(arch, imap, "Failed to read index map in formic_ExpandAPJ_command_function_detail.");

  // get small jastrow factor data
  const std::vector<S> & sm_aa_vec = wfn.get_wf(sm_aa_name)->data_vec();
  const std::vector<S> & sm_bb_vec = wfn.get_wf(sm_bb_name)->data_vec();
  const std::vector<S> & sm_ab_vec = wfn.get_wf(sm_ab_name)->data_vec();

  // get large jastrow factor data
  std::vector<S> & lg_aa_vec = wfn.get_wf(lg_aa_name)->data_vec();
  std::vector<S> & lg_bb_vec = wfn.get_wf(lg_bb_name)->data_vec();
  std::vector<S> & lg_ab_vec = wfn.get_wf(lg_ab_name)->data_vec();

  // make sure large set is at least as large as small set
  if ( sm_norb < 2       ) throw formic::Exception("sm_norb must be >= 2 in formic_ExpandAPJ_command_function_detail");
  if ( lg_norb < sm_norb ) throw formic::Exception("lg_norb must be >= sm_norb in formic_ExpandAPJ_command_function_detail");

  // check sizes
  if ( sm_aa_vec.size() != 4*sm_norb*(sm_norb-1)/2 ) throw formic::Exception("sm_aa_vec is the wrong size in formic_ExpandAPJ_command_function_detail");
  if ( sm_bb_vec.size() != 4*sm_norb*(sm_norb-1)/2 ) throw formic::Exception("sm_bb_vec is the wrong size in formic_ExpandAPJ_command_function_detail");
  if ( sm_ab_vec.size() != 4*sm_norb*sm_norb       ) throw formic::Exception("sm_ab_vec is the wrong size in formic_ExpandAPJ_command_function_detail");
  if ( lg_aa_vec.size() != 4*lg_norb*(lg_norb-1)/2 ) throw formic::Exception("lg_aa_vec is the wrong size in formic_ExpandAPJ_command_function_detail");
  if ( lg_bb_vec.size() != 4*lg_norb*(lg_norb-1)/2 ) throw formic::Exception("lg_bb_vec is the wrong size in formic_ExpandAPJ_command_function_detail");
  if ( lg_ab_vec.size() != 4*lg_norb*lg_norb       ) throw formic::Exception("lg_ab_vec is the wrong size in formic_ExpandAPJ_command_function_detail");

  // check index map size
  if ( imap.size() != sm_norb )
    throw formic::Exception("size of index map is %i, but should be %i in formic_ExpandAPJ_command_function_detail") % imap.size() % sm_norb;

  // check index map contents
  for (std::vector<int>::const_iterator it = imap.begin(); it != imap.end(); it++) {
    if ( *it < 0 || *it >= lg_norb )
      throw formic::Exception("out of range value \"%i\" in index map in formic_ExpandAPJ_command_function_detail") % *it;
    for (std::vector<int>::const_iterator jt = imap.begin(); jt != it; jt++)
      if ( *jt == *it )
        throw formic::Exception("repeated index \"%i\" in index map in formic_ExpandAPJ_command_function_detail") % *jt;
  }

  // get ioff array
  std::vector<int> ioff;
  formic::get_pair_ioff(lg_norb, ioff);

  // expand aa
  for (int j = 1; j < sm_norb; j++)
  for (int i = 0; i < j; i++) {
    const int p = imap[i];
    const int q = imap[j];
    if ( p < q ) {
      const int sm_pos = 4 * ( i + ioff[j] );
      const int lg_pos = 4 * ( p + ioff[q] );
      lg_aa_vec[lg_pos+0] += sm_aa_vec[sm_pos+0];
      lg_aa_vec[lg_pos+1] += sm_aa_vec[sm_pos+1];
      lg_aa_vec[lg_pos+2] += sm_aa_vec[sm_pos+2];
      lg_aa_vec[lg_pos+3] += sm_aa_vec[sm_pos+3];
    } else {
      const int sm_pos = 4 * ( i + ioff[j] );
      const int lg_pos = 4 * ( q + ioff[p] );
      lg_aa_vec[lg_pos+0] += sm_aa_vec[sm_pos+0];
      lg_aa_vec[lg_pos+2] += sm_aa_vec[sm_pos+1];
      lg_aa_vec[lg_pos+1] += sm_aa_vec[sm_pos+2];
      lg_aa_vec[lg_pos+3] += sm_aa_vec[sm_pos+3];
    }
  }

  // expand bb
  for (int j = 1; j < sm_norb; j++)
  for (int i = 0; i < j; i++) {
    const int p = imap[i];
    const int q = imap[j];
    if ( p < q ) {
      const int sm_pos = 4 * ( i + ioff[j] );
      const int lg_pos = 4 * ( p + ioff[q] );
      lg_bb_vec[lg_pos+0] += sm_bb_vec[sm_pos+0];
      lg_bb_vec[lg_pos+1] += sm_bb_vec[sm_pos+1];
      lg_bb_vec[lg_pos+2] += sm_bb_vec[sm_pos+2];
      lg_bb_vec[lg_pos+3] += sm_bb_vec[sm_pos+3];
    } else {
      const int sm_pos = 4 * ( i + ioff[j] );
      const int lg_pos = 4 * ( q + ioff[p] );
      lg_bb_vec[lg_pos+0] += sm_bb_vec[sm_pos+0];
      lg_bb_vec[lg_pos+2] += sm_bb_vec[sm_pos+1];
      lg_bb_vec[lg_pos+1] += sm_bb_vec[sm_pos+2];
      lg_bb_vec[lg_pos+3] += sm_bb_vec[sm_pos+3];
    }
  }

  // expand ab
  for (int j = 0; j < sm_norb; j++)
  for (int i = 0; i < sm_norb; i++) {
    const int p = imap[i];
    const int q = imap[j];
    const int sm_pos = 4 * ( i + j * sm_norb );
    const int lg_pos = 4 * ( p + q * lg_norb );
    lg_ab_vec[lg_pos+0] += sm_ab_vec[sm_pos+0];
    lg_ab_vec[lg_pos+1] += sm_ab_vec[sm_pos+1];
    lg_ab_vec[lg_pos+2] += sm_ab_vec[sm_pos+2];
    lg_ab_vec[lg_pos+3] += sm_ab_vec[sm_pos+3];
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Converts a set of non-simple jastrows into a set of simple jastrows.
///
/// \param[in]       userinp  the user's input options
/// \param[in,out]   wfn      the wavefunction object holding the jastrow factor data
/// \param[in,out]   arch     archive to read command arguments from
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S>
static void formic_ToSimpleJastrows_command_function_detail(const formic::InputBase & userinp,
                                                            formic::Wavefunction<S> & wfn,
                                                            formic::Archive & arch) {

  throw formic::Exception("ToSimpleJastrows command is not working quite right.  You should fix it!");

  // variables to describe the two sets of jastrow factors
  std::string n_aa_name, n_bb_name, n_ab_name, s_aa_name, s_bb_name, s_ab_name;

  // read in the non-simple jastrow factor names
  formic::mpi::read_and_bcast(arch, n_aa_name, "Failed to read n_aa_name in formic_ToSimpleJastrows_command_function_detail.");
  formic::mpi::read_and_bcast(arch, n_bb_name, "Failed to read n_bb_name in formic_ToSimpleJastrows_command_function_detail.");
  formic::mpi::read_and_bcast(arch, n_ab_name, "Failed to read n_ab_name in formic_ToSimpleJastrows_command_function_detail.");

  // read in the simple jastrow factor names
  formic::mpi::read_and_bcast(arch, s_aa_name, "Failed to read s_aa_name in formic_ToSimpleJastrows_command_function_detail.");
  formic::mpi::read_and_bcast(arch, s_bb_name, "Failed to read s_bb_name in formic_ToSimpleJastrows_command_function_detail.");
  formic::mpi::read_and_bcast(arch, s_ab_name, "Failed to read s_ab_name in formic_ToSimpleJastrows_command_function_detail.");

  // get simple jastrow factor data
  std::vector<S> & s_aa_vec = wfn.get_wf(s_aa_name)->data_vec();
  std::vector<S> & s_bb_vec = wfn.get_wf(s_bb_name)->data_vec();
  std::vector<S> & s_ab_vec = wfn.get_wf(s_ab_name)->data_vec();

  // get non-simple jastrow factor data
  const std::vector<S> & n_aa_vec = wfn.get_wf(n_aa_name)->data_vec();
  const std::vector<S> & n_bb_vec = wfn.get_wf(n_bb_name)->data_vec();
  const std::vector<S> & n_ab_vec = wfn.get_wf(n_ab_name)->data_vec();

  // zero the simple jastrow data
  std::fill(s_aa_vec.begin(), s_aa_vec.end(), formic::zero(S()));
  std::fill(s_bb_vec.begin(), s_bb_vec.end(), formic::zero(S()));
  std::fill(s_ab_vec.begin(), s_ab_vec.end(), formic::zero(S()));

  // get number of orbitals
  const int norb = userinp.get<int>("nao");

  // do conversion
  if ( userinp.get<bool>("restrict_jastrow") ) {

    formic::Exception("formic_ToSimpleJastrows_command_function_detail not yet implemented for restricted jastrows");

  } else {

    // check sizes
    if ( n_aa_vec.size() != 4*norb*(norb-1)/2 ) throw formic::Exception("n_aa_vec is the wrong size in formic_ToSimpleJastrows_command_function_detail");
    if ( n_bb_vec.size() != 4*norb*(norb-1)/2 ) throw formic::Exception("n_bb_vec is the wrong size in formic_ToSimpleJastrows_command_function_detail");
    if ( n_ab_vec.size() != 4*norb*norb       ) throw formic::Exception("n_ab_vec is the wrong size in formic_ToSimpleJastrows_command_function_detail");
    if ( s_aa_vec.size() !=   norb*(norb+1)/2 ) throw formic::Exception("s_aa_vec is the wrong size in formic_ToSimpleJastrows_command_function_detail");
    if ( s_bb_vec.size() !=   norb*(norb+1)/2 ) throw formic::Exception("s_bb_vec is the wrong size in formic_ToSimpleJastrows_command_function_detail");
    if ( s_ab_vec.size() !=   norb*norb       ) throw formic::Exception("s_ab_vec is the wrong size in formic_ToSimpleJastrows_command_function_detail");

    // get ioff array
    std::vector<int> ioff;
    formic::get_pair_ioff(norb, ioff);

    // expand aa
    for (int j = 1; j < norb; j++)
    for (int i = 0; i < j; i++) {

      const int n_pos = 4 * ( i + ioff.at(j) );

      s_aa_vec.at( formic::cmpd_pair_index_incl(i, i, norb) ) -= n_aa_vec.at(n_pos+0); // ni
      s_aa_vec.at( formic::cmpd_pair_index_incl(j, j, norb) ) -= n_aa_vec.at(n_pos+0); // nj
      s_aa_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) += n_aa_vec.at(n_pos+0); // ni nj

      s_aa_vec.at( formic::cmpd_pair_index_incl(j, j, norb) ) += n_aa_vec.at(n_pos+1); // nj
      s_aa_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) -= n_aa_vec.at(n_pos+1); // ni nj

      s_aa_vec.at( formic::cmpd_pair_index_incl(i, i, norb) ) += n_aa_vec.at(n_pos+2); // ni
      s_aa_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) -= n_aa_vec.at(n_pos+2); // ni nj

      s_aa_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) += n_aa_vec.at(n_pos+3); // ni nj

    }

    // expand bb
    for (int j = 1; j < norb; j++)
    for (int i = 0; i < j; i++) {

      const int n_pos = 4 * ( i + ioff.at(j) );

      s_bb_vec.at( formic::cmpd_pair_index_incl(i, i, norb) ) -= n_bb_vec.at(n_pos+0); // ni
      s_bb_vec.at( formic::cmpd_pair_index_incl(j, j, norb) ) -= n_bb_vec.at(n_pos+0); // nj
      s_bb_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) += n_bb_vec.at(n_pos+0); // ni nj

      s_bb_vec.at( formic::cmpd_pair_index_incl(j, j, norb) ) += n_bb_vec.at(n_pos+1); // nj
      s_bb_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) -= n_bb_vec.at(n_pos+1); // ni nj

      s_bb_vec.at( formic::cmpd_pair_index_incl(i, i, norb) ) += n_bb_vec.at(n_pos+2); // ni
      s_bb_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) -= n_bb_vec.at(n_pos+2); // ni nj

      s_bb_vec.at( formic::cmpd_pair_index_incl(i, j, norb) ) += n_bb_vec.at(n_pos+3); // ni nj

    }

    // expand ab
    for (int j = 0; j < norb; j++)
    for (int i = 0; i < norb; i++) {

      const int n_pos = 4 * ( i + j * norb );

      s_aa_vec.at( formic::cmpd_pair_index_incl(i, i, norb) ) -= n_ab_vec.at(n_pos+0); // ni
      s_bb_vec.at( formic::cmpd_pair_index_incl(j, j, norb) ) -= n_ab_vec.at(n_pos+0); // nj
      s_ab_vec.at( i * norb + j                             ) += n_ab_vec.at(n_pos+0); // ni nj

      s_bb_vec.at( formic::cmpd_pair_index_incl(j, j, norb) ) += n_ab_vec.at(n_pos+1); // nj
      s_ab_vec.at( i * norb + j                             ) -= n_ab_vec.at(n_pos+1); // ni nj

      s_aa_vec.at( formic::cmpd_pair_index_incl(i, i, norb) ) += n_ab_vec.at(n_pos+2); // ni
      s_ab_vec.at( i * norb + j                             ) -= n_ab_vec.at(n_pos+2); // ni nj

      s_ab_vec.at( i * norb + j                             ) += n_ab_vec.at(n_pos+3); // ni nj

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "ExpandAPJ" command, which adds a set of all pair jastrow
///          factors from a smaller set of orbitals to a larger set of orbitals containing
///          the original smaller set.
///
///          The first four arguments are the number of smaller set orbitals and the names of the
///          smaller aa, bb, and ab jastrow factor data objects.
///
///          The next four arguments are the number of larger set orbitals and the names of the
///          larger aa, bb, and ab jastrow factor data objects.
///
///          The last argument is a vector giving the position of each smaller-orbital-set orbital
///          in the ordering of the larger orbital set.
///
///          For example, expanding the variables saa, sbb, and sab of a 2-orbital all pair
///          jastrow into the variables laa, lbb, and lab of an 8-orbital all pair jastrow when
///          the two small-orbital-set orbitals 0 and 1 map to orbitals 3 and 5 would use the
///          command:
///
///            ExpandAPJ  2 saa sbb sab  8 laa lbb lab  [ 3 5 ]
///
/// \param[in,out]   arch     archive to read command arguments from
/// \param[in,out]   data     data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_ExpandAPJ_command_function(formic::Archive & arch, formic::AllData & data) {

  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_ExpandAPJ_command_function_detail(data.userinp, data.wd, arch);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_ExpandAPJ_command_function_detail(data.userinp, data.wc, arch);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" when processing ExpandAPJ command.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "ToSimpleJastrows" command, which converts a set of non-simple
///          jastrow variables into simple jastrow format.
///
///          The first three arguments are the names of the non-simple jastrow data objects.
///
///          The second three arguments are the names of the data objects to place the simple
///          jastrows in.
///
///          Example command:
///
///            ToSimpleJastrows  jaa jbb jab  saa sbb sab
///
/// \param[in,out]   arch     archive to read command arguments from
/// \param[in,out]   data     data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_ToSimpleJastrows_command_function(formic::Archive & arch, formic::AllData & data) {

  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_ToSimpleJastrows_command_function_detail(data.userinp, data.wd, arch);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_ToSimpleJastrows_command_function_detail(data.userinp, data.wc, arch);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" when processing ToSimpleJastrows command.") % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Adds a partial particle number projection to the specified jastrow factors.
///
/// \param[in]       userinp  the user's input options
/// \param[in,out]   wfn      the wavefunction object holding the jastrow factor data
/// \param[in,out]   arch     archive to read command arguments from
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S>
static void formic_JastNumProj_command_function_detail(const formic::InputBase & userinp,
                                                       formic::Wavefunction<S> & wfn,
                                                       formic::Archive & arch) {

  // get whether or not jastrow factor is restricted
  const bool rest_jast = userinp.get<bool>("restrict_jastrow");

  // the spin of the electrons to project (0 == alpha, 1 == beta, 2 == both)
  int proj_spin;
  formic::mpi::read_and_bcast(arch, proj_spin, "Failed to read proj_spin in formic_JastNumProj_command_function_detail.");
  if ( formic::mpi::rank() == 0 && ( proj_spin < 0 || proj_spin > 2 ) )
    throw formic::Exception("first argument to JastNumProj is the electron spin and it must be either 0, 1, or 2");

  // the projection strength
  double alpha = 0.0;
  formic::mpi::read_and_bcast(arch, alpha, "Failed to read alpha in formic_JastNumProj_command_function_detail.");
  if ( formic::mpi::rank() == 0 && alpha < 0.0 )
    throw formic::Exception("second argument to JastNumProj is the projection strength and it must be >= 0.0");

  // the desired number of electrons to project to
  int nproj;
  formic::mpi::read_and_bcast(arch, nproj, "Failed to read nproj in formic_JastNumProj_command_function_detail.");
  if ( formic::mpi::rank() == 0 && nproj < 0 )
    throw formic::Exception("third argument to JastNumProj is the number of particles to project to and it must be >= 0");

  // read in the jastrow factor names
  std::string aa_name, bb_name, ab_name;
  formic::mpi::read_and_bcast(arch, aa_name, "Failed to read aa_name (4th arg) in formic_JastNumProj_command_function_detail.");
  if ( !rest_jast )
    formic::mpi::read_and_bcast(arch, bb_name, "Failed to read bb_name (5th arg) in formic_JastNumProj_command_function_detail.");
  formic::mpi::read_and_bcast(arch, ab_name, "Failed to read ab_name (6th arg) in formic_JastNumProj_command_function_detail.");

  // read in the set of orbitals on which to project
  std::vector<int> proj_orb;
  formic::mpi::read_and_bcast(arch, proj_orb, "Failed to read proj_orb (7th arg) in formic_JastNumProj_command_function_detail.");

  // get jastrow factor data
  std::vector<S> & aa_vec = wfn.get_wf(aa_name)->data_vec();
  std::vector<S> & bb_vec = ( rest_jast ? wfn.get_wf(aa_name)->data_vec() : wfn.get_wf(bb_name)->data_vec() );
  std::vector<S> & ab_vec = wfn.get_wf(ab_name)->data_vec();

  // check sizes
  const int norb = userinp.get<int>("nao");
  if ( userinp.get<bool>("simple_jastrow") ) {
    const int corr_ab_size = ( rest_jast ? norb*(norb+1)/2 : norb*norb );
    if ( aa_vec.size() != norb*(norb+1)/2 ) throw formic::Exception("aa_vec is the wrong size in formic_JastNumProj_command_function_detail");
    if ( !rest_jast )
      if ( bb_vec.size() != norb*(norb+1)/2 ) throw formic::Exception("bb_vec is the wrong size in formic_JastNumProj_command_function_detail");
    if ( ab_vec.size() != corr_ab_size    ) throw formic::Exception("ab_vec is the wrong size in formic_JastNumProj_command_function_detail");
  } else {
    if ( aa_vec.size() != 4*norb*(norb-1)/2 ) throw formic::Exception("aa_vec is the wrong size in formic_JastNumProj_command_function_detail");
    if ( !rest_jast )
      if ( bb_vec.size() != 4*norb*(norb-1)/2 ) throw formic::Exception("bb_vec is the wrong size in formic_JastNumProj_command_function_detail");
    if ( ab_vec.size() != 4*norb*norb       ) throw formic::Exception("ab_vec is the wrong size in formic_JastNumProj_command_function_detail");
  }

  // check projection orbitals
  if ( formic::mpi::rank() == 0 && proj_spin == 2 && ( 2*proj_orb.size() < std::max(1,nproj) || proj_orb.size() > norb ) )
    throw formic::Exception("invalid length \"%i\" of vector of projection orbitals in formic_JastNumProj_command_function_detail") % proj_orb.size();
  if ( formic::mpi::rank() == 0 && proj_spin  < 2 && ( proj_orb.size() < std::max(1,nproj) || proj_orb.size() > norb ) )
    throw formic::Exception("invalid length \"%i\" of vector of projection orbitals in formic_JastNumProj_command_function_detail") % proj_orb.size();
  for (int i = 0; i < proj_orb.size(); i++)
    if ( formic::mpi::rank() == 0 && ( proj_orb.at(i) < 0 || proj_orb.at(i) >= norb ) )
      throw formic::Exception("bad value \"%i\" for projection orbital %i in formic_JastNumProj_command_function_detail") % proj_orb.at(i) % i;

  // sort projection orbitals
  std::sort(proj_orb.begin(), proj_orb.end());

  // get ioff array
  std::vector<int> ioff;
  formic::get_pair_ioff(norb, ioff);

  // ------ Projection for total electron number (alpha and beta) in given orbitals ------
  if ( proj_spin == 2 ) {

    // add the projection terms to the jastrow factors
    if ( userinp.get<bool>("simple_jastrow") ) {

      // get numbers of alpha, beta, and total electrons
      const int na = userinp.get<int>("nalpha");
      const int nb = userinp.get<int>("nbeta");
      const int nt = na + nb;

      // constant term
      for (int p = 0; p < norb; p++)
        aa_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) -= formic::unity(S()) * ( alpha * nproj * nproj / double(nt) );
      if ( ! rest_jast )
        for (int p = 0; p < norb; p++)
          bb_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) -= formic::unity(S()) * ( alpha * nproj * nproj / double(nt) );

      // loop over first projected orbital index
      for (int j = 0; j < proj_orb.size(); j++) {

        // get first projected orbital index
        const int p = proj_orb[j];

        // cross term
        aa_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) += formic::unity(S()) * ( 2.0 * alpha * nproj );
        if ( ! rest_jast )
          bb_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) += formic::unity(S()) * ( 2.0 * alpha * nproj );

        // loop over second projected orbital index
        for (int i = 0; i < proj_orb.size(); i++) {

          // get second projected orbital index
          const int q = proj_orb[i];

          // same-spin double-number-operator terms
          aa_vec.at( formic::cmpd_pair_index_incl(std::min(p,q), std::max(p,q), norb) ) -= formic::unity(S()) * alpha;
          if ( ! rest_jast )
            bb_vec.at( formic::cmpd_pair_index_incl(std::min(p,q), std::max(p,q), norb) ) -= formic::unity(S()) * alpha;

          // opposite-spin double-number-operator terms
          ab_vec.at( formic::fqmc::jastrow_ab_cmpd_ind(rest_jast, p, q, norb) )
                           -= formic::unity(S()) * alpha * ( rest_jast && p != q ? 1.0 : 2.0 );

        }

      }

    } else {

      throw formic::Exception("projection in total electron number (alpha and beta) only implemented for simple jastrows");

    }

  // ------ Projection for spin-specific electron number (alpha or beta) in given orbitals ------
  } else {

    // get aa or bb jastrow factor data array
    std::vector<S> & xx_vec = ( proj_spin == 0 ? aa_vec : bb_vec );

    // add the projection terms to the jastrow factors
    if ( userinp.get<bool>("simple_jastrow") ) {

      // get number of alpha or beta electrons
      const int nss = ( proj_spin == 0 ? userinp.get<int>("nalpha") : userinp.get<int>("nbeta") );

      // constant term
      for (int p = 0; p < norb; p++)
        xx_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) -= formic::unity(S()) * ( alpha * nproj * nproj / double(nss) );

      // loop over first projected orbital index
      for (int j = 0; j < proj_orb.size(); j++) {

        // get first projected orbital index
        int p = proj_orb[j];

        // same-number double-number-operator term
        xx_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) -= formic::unity(S()) * alpha;

        // cross term
        xx_vec.at( formic::cmpd_pair_index_incl(p, p, norb) ) += formic::unity(S()) * ( 2.0 * alpha * nproj );

        // loop over second projected orbital index
        for (int i = 0; i < j; i++) {

          // get second projected orbital index
          int q = proj_orb[i];

          // different-number double-number-operator term
          xx_vec.at( formic::cmpd_pair_index_incl(q, p, norb) ) -= formic::unity(S()) * 2.0 * alpha;

        }

      }

    } else {

      //const S c10 = formic::unity(S()) * 0.5 * double(2*nproj-1) * alpha / double(proj_orb.size()-1);
      //const S c10 = formic::unity(S()) * 2.0 * alpha * double(nproj-1) / double(proj_orb.size()-1);
      const S c10 = formic::unity(S()) * alpha * double(2*nproj-1) / double(proj_orb.size()-1);
      const S c11 = 2.0 * c10 - 2.0 * alpha * formic::unity(S());
      for (int j = 1; j < proj_orb.size(); j++)
      for (int i = 0; i < j; i++) {

        // get the location of the specific jastrow factor in the data array
        int p = proj_orb[i];
        int q = proj_orb[j];
        if ( p > q )
          std::swap(p,q);
        const int pos = 4 * ( p + ioff[q] );

        xx_vec[pos+1] += c10;
        xx_vec[pos+2] += c10;
        xx_vec[pos+3] += c11;

      }

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Function called for the "JastNumProj" command, which adds to the specified jastrow
///          factors a partial number projection of electrons of the given spin on the given set
///          of orbitals.
///          The projection is partial in the sense that it is only a strict projection in the
///          limit as the parameter alpha heads to positive infinity.
///
///          The first argument is the spin of the electrons to operate on, 0 for alpha,
///          1 for beta, or 2 for both.
///
///          The second argument is the strength of the partial projection, given by alpha,
///          which more positive alphas giving more strict projections.
///
///          The third argument is the desired number of electrons to project towards.
///
///          The 4th, 5th, and (possibly) 6th arguments are the names of the aa, bb, and ab
///          all pair jastrow factor data objects to which to add the projection.
///          For restricted jastrow factors, the bb argument should be omitted.
///
///          The final argument is a vector of orbitals on which to perform the projection.
///
///          For example, projecting towards 4 beta electrons on orbitals 0, 1, 3, 4, 5, and 6
///          with a projection strength of 2.5 and jastrow factors named jaa, jbb, and jab
///          would use the command:
///
///            JastNumProj  1  2.5  4  jaa  jbb  jab  [ 0 1 3 4 5 6 ]
///
///          For a restricted jastrow example, projecting towards 3 total electrons on orbitals
///          0, 1, and 4 with a projection strength of 0.5 and jastrow factors named jaa and jab
///          would use the command:
///
///            JastNumProj  2  0.5  3  jaa  jab  [ 0 1 4 ]
///
/// \param[in,out]   arch     archive to read command arguments from
/// \param[in,out]   data     data on which to execute the command
///
///////////////////////////////////////////////////////////////////////////////////////////////////
static void formic_JastNumProj_command_function(formic::Archive & arch, formic::AllData & data) {

  if (data.userinp.get<std::string>("scalar_type") == "real") {
    formic_JastNumProj_command_function_detail(data.userinp, data.wd, arch);
  } else if (data.userinp.get<std::string>("scalar_type") == "complex") {
    formic_JastNumProj_command_function_detail(data.userinp, data.wc, arch);
  } else 
    throw formic::Exception("Unknown scalar type \"%s\" when processing JastNumProj command.")
            % data.userinp.get<std::string>("scalar_type");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Add restricted geminal power related commands to the formic command system
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::init_all_pair_jast_commands() {

  // add command to construct an all pair jastrow coefficient object
  formic::add_command("AllPairJast", &formic_AllPairJast_command_function);

  // add command to expand a set of all pair jastrows into a larger orbital space
  formic::add_command("ExpandAPJ", &formic_ExpandAPJ_command_function);

  // add command to add a partial number projection to a set of jastrow factors
  formic::add_command("JastNumProj", &formic_JastNumProj_command_function);

  // add command to convert non-simple to simple jastrows
  formic::add_command("ToSimpleJastrows", &formic_ToSimpleJastrows_command_function);

}
